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We study the quantum version of a lattice model whose classical counterpart captures the physics 
of structural glasses. We discuss the role of quantum fluctuations in such systems and in particular 
their interplay with the amorphous order developed in the glass phase. We show that quantum 
fluctuations might facilitate the formation of the glass at low enough temperature. We also show 
that the glass transition becomes a first-order transition between a superfluid and an insulating glass 
at very low temperature, and is therefore accompanied by phase coexistence between superfluid and 
glassy regions. 

I. INTRODUCTION 

It is a well established fact that certain liquids form a glass if cooled sufficiently fast to avoid crystallization. 
The simplest examples of such glass- formers are models of particles interacting via Lennard- Jones-like [1, 2] and hard 
spheres [3] potentials, but a similar phenomenology is observed in many real systems like Silica, molecular systems such 
as Glycerol and OTP, and many polymers [4, 5] . The precise meaning of "sufficiently fast to avoid crystallization" is 
strongly system-dependent, with time scales ranging from nanoseconds for three-dimensional monoatomic systems, to 
hours or centuries for more complex systems like binary mixtures or molecular and polymeric systems. If crystallization 
is avoided, below the melting temperature T m the system enters a "supercooled liquid" state [4] , which can be treated as 
a metastable equilibrium state and hence described using tools from equilibrium statistical mechanics. The transition 
to the glass happens at a lower temperature T g , at which the relaxation time of the metastable liquid exceeds the 
experimentally accessible time scales and the system behaves as a solid for any practical purpose [4, 5]. The glass 
transition can be driven by lowering the temperature T or increasing the density p, and one can define a line of 
transitions T g (p) which delimits the "glassy" region of the phase diagram. It is interesting to note that although the 
glass transition happens out of equilibrium, its phenomenology is very close to that of a thermodynamical second 
order phase transition, signaled by a jump in specific heat in absence of any latent heat of solidification. It has been 
conjectured by many authors that a true thermodynamic phase transition happens at a lower temperature Tx(p) and 
is responsible for the non-equilibrium glass transition at T g . 

It is reasonable to expect that at low enough temperatures and high enough density (the precise scales depending on 
the details of the interaction potential), quantum fluctuations (quantified for instance by the ratio T = A/a, where A 
is the De Broglie thermal wavelength and a is the interparticle distance) becomes relevant. A typical example is that 
of He 4 , that forms a solid at temperatures so low that quantum effects are very important [6]. Another interesting 
example is that of hydrogen at very high pressure and low temperature [7]. Yet, the possible presence of glassy 
phases in these quantum systems has been proposed only very recently [8, 9], and its experimental observation is still 
debated [10, 11]. These preliminary observations raise the natural question whether the classical picture of the glass 
transition is strongly modified by quantum fluctuations. In particular, one would like to know what is the nature 
of the transition (does it remain second order?), what is the shape of the glass transition line in the (T,p,T) space, 
what happens to the quantum dynamics in the vicinity of the glass transition (is tunnelling between distinct glassy 
metastable states possible?), and so on. Due to the very complex nature of the glass transition, that is characterized 
by the breaking of translational invariance towards an amorphous ordering and by very slow dynamics, it is very hard 
to address these questions using numerical methods such as Quantum Monte Carlo [8] or exact diagonalization. 

As already mentioned, for classical systems a widespread approach to the glass transition consists in treating it as 
a true phase transition, using tools of equilibrium statistical mechanics. After several decades of efforts, a coherent 
theory of the glass transition based on this approach has emerged and now goes under the name of Random First 
Order Transition (RFOT) theory [12]. According to this theory, the glass transition belongs to a new classification: it 
is a second order transition from a thermodynamic point of view, but its order parameter (that we will define later on) 
jumps at the transition and many physical properties around the transition are related to some nucleation-like physics 
(hence the transition looks more like a first order transition in this respect, which is at the origin of the name RFOT). 
The main inspiration for this theory came from density functional theory [13, 14], from Mode-Coupling Theory (a 
dynamic theory of supercooled liquids that emerged in the context of critical phenomena) [14, 15] and from the exact 
solution of a class of mean-field spin glass models (the so-called p-spin glasses) [16-18] that display a phenomenology 
extremely close to that of glass-forming liquids. A complete review of RFOT theory goes well beyond the scope of 
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this paper, and we refer the reader to the excellent reviews [12, 19, 20] for more details. What is relevant for the 
present discussion is the existence of mean field lattice models of glasses [21-23] that are extremely simple (they are 
formulated in terms of hard particles on a lattice with simple local many body repulsion), exactly solvable (thanks 
to the mean field nature of the lattice one consider) and at the same time display the full RFOT phenomenology: in 
particular for these models one can compute the glass transition line Tx(p), and one can also study the dynamics on 
approaching the glass transition from the liquid. These models are a very natural starting point for the investigation 
of the glass transition in the quantum case. Moreover, the RFOT construction can be implemented under controlled 
approximations to compute quantitatively the glass transition in realistic systems, such as Lennard- Jones [24, 25] or 
hard spheres [26]. 



The extension of RFOT theory to quantum systems has been initiated long time ago by studying the quantum 
version of the p-spin glasses [27-31], stripe glasses [32] and of some optimization problems [33] that share the same 
qualitative phenomenology. The main result of these studies is that the high temperature (RFOT- like) glass transition 
is turned into a first order transition between the liquid and the glass at very low temperature. However, the picture 
provided by these studies is not complete, since all the models that have been studied up to now are in the glass phase 
at low enough T in absence of quantum fluctuations; in other words, none of these models displays, in the classical 
limit r = 0, a glass transition at low temperature as a function of another control parameter (e.g. the density). The 
possibility of studying the glass transition at low T as a function of p adds an important ingredient to the discussion, 
which is entropy. In fact, the models studied in [27-33] have a non-extensively degenerate ground state, hence their 
entropy at T = vanishes. Conversely, lattice glass models might have a finite entropy even at T = (think for 
instance to hard spheres [26]) and in this case the glass transition is completely driven by entropy. A preliminary 
attempt to take into account the role of entropy in this class of problems has been reported for a toy model in [34] , 
where it has been shown that due to the interplay of classical entropy and quantum fluctuations, the latter may induce 
the formation of a glass: from a liquid phase at V = the system can fall in a glassy situation as soon as T > 0. A 
similar phenomenon has been found in [35] for a system of quantum hard spheres. 

Finally, a very interesting point that has been raised in the study of quantum glasses is the possibility that a glass 
display superfluidity [8, 9]. This would lead to a particular supersolid phase (a "superglass"). However, the existence 
of such a phase has by now only been shown numerically [8] or analytically in a quite special model [9]. For systems 
displaying a standard second order spin glass transition [36] (which is quite different from a RFOT), a superglass 
phase has been shown to exist generically [37, 38]. However, the physics of RFOT is quite different from that of a 
second order spin glass transition, therefore the extension of these results to the RFOT scenario is not straightforward. 



In this work we shall study a Bosonic quantum extension of the Biroli-Mezard lattice glass model [21]. Our aim 
is to study the effect of quantum fluctuations in this system and in particular their interplay with the amorphous 
order developed in the glass phase, and with the large degeneracy of states leading to a finite entropy at T = 0. 
Specifically, we want to compute the full phase diagram of the system to find the liquid-glass transition line, and look 
for a superglass phase. We will do this by means of the recently developed quantum cavity method [39-41]. We find 
a complex phase diagram at low temperature, characterized by two unexpected phenomena: i) A re-entrance of the 
glass transition line towards lower density on increasing quantum fluctuations, driven by entropy; ii) A first order 
phase transition between superfluid and glass at zero temperature. Our results suggest that a true superglass phase is 
not present in the model, but we can establish the presence of phase coexistence between superfluid and glass, which 
might have interesting phenomenological consequences. 

The rest of the paper is divided as follows. In section II we introduce the classical model clarifying its importance in 
the description of classical glasses. In section II A we derive the so-called "cavity-equations" describing the system in 
the replica symmetric ansatz, which provides the exact solution of the liquid phase. Then, in section II B and II C, we 
describe how to access its glass phase within the replica symmetry breaking scheme. A summary of the results for the 
classical model is given in section II D. In section III we generalize the model as a quantum many-body system, made 
of hard-core bosons, with the same interaction potential as the classical one but where quantum effects are induced 
by an hopping term. The definition of the model is given in section III A. The rest of Sec. Ill is more technical and 
can be skipped by readers not interested in the formalism of the quantum cavity method used to solve the model; 
in section III B we derive the "cavity equations" which describe the physics of the system and in sections III C and 
HID we sketch their numerical resolution. In section IV we present the main results for the quantum system. In 
section IV A we define the order parameters of the different phases; in section IV B we discuss the details of the cavity 
computations; in section IV C we report the main results on the phase diagram of the model. In section IV D we 
investigate the quantum dynamics of the model: we show that at the glass transition there is ergodicity breaking 
and we discuss its relation with the spectrum of the quantum Hamiltonian. In section IV E we discuss some exact 
diagonalization results. Finally, in section V we summarize and discuss our results. 
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II. THE CLASSICAL MODEL 

The classical model considered in the following is a generalization of a lattice glass model proposed in [21, 22] and 
widely used as a prototype for the description of particle systems undergoing a glass transition. Its degrees of freedom 
are occupation numbers n, = {0,1} of the sites i = 1, . . . , N of a graph, a particle being present on site i if and only 
if rii = 1. If a site i is occupied it feels a repulsive interaction with the particles on the neighboring sites of the graph, 
the set of which shall be denoted di. More precisely, each particle can have at most £ neighbors without any energy 
cost and it is subject to an interaction V > for every neighboring particle in excess. The Hamiltonian describing 
this interaction is the following: 

H = V^rn qi 6(qi) - p^rn , (1) 

i i 

where qi = J2jedi n i ~ ^ ^ s tnc numDCr of neighbors "in excess". Here and in the following 9{x) — if x < and 
9(x) = 1 if x > 0. We included in H the chemical potential p, conjugated to the total particle number. The original 
model of [21, 22] corresponds to the case V — > oo; configurations where one particle has more than I neighboring sites 
occupied are then strictly forbidden. This finite V version will be more convenient for the quantum generalization 
described later on. 

The model above can be defined on any finite graph. Its finite dimensional version has been studied numerically by 
Monte Carlo simulations in [21, 42, 43], demonstrating that it exhibits the phenomenology of glassy systems. Note 
that there is no quenched disorder in the Hamiltonian; as in real liquids the disorder and frustration that yield an 
amorphous glassy order are in fact self-generated by the different possible arrangements of particles. In order to 
obtain analytical results we will study the model on the Bcthe lattice, as in [21, 22]. This underlying geometry is 
known to give a mean field description of the system, with the advantage of preserving the notion of distance and that 
of finite connectivity [41, 44]. More precisely, Bethe lattices are intended here as random c-regular graphs, i.e. graphs 
chosen uniformly among the set of graphs with N vertices where all sites have precisely c neighbors. These graphs 
converge locally to trees in the thermodynamic (large N) limit. There are several motivations to study the problem 
on this geometry [44]. Due to the local tree-like structure, in fact, on this kind of graphs one is able to take explicitly 
the thermodynamic limit and derive the exact solution of the model. On the other hand, having no boundary and 
all sites playing statistically the same role, random regular graphs allow to get rid of boundary problems which affect 
trees. Still, even if locally tree-like, they do have loops, which are an essential ingredient for the frustration induced 
in the system. Moreover, their randomness forbids crystalline ordering, making easier the study of glassy phases. All 
these properties make the Bethe lattice a convenient playground for the study of frustrated systems. More generally, 
however, the Bcthe lattice can also give insight into Euclidean lattices. In fact, the solution that one gets on random 
regular graphs can be thought of as a "refined" (with respect to the standard fully-connected model) mean field 
approximation for a finite dimensional lattice with the same connectivity. Finally the free energy of the Bethe lattice 
has the important property of being self-averaging, in the thermodynamic limit, with respect to the uniform measure 
which defines the ensemble. This ensures that for large enough size already one exemplar is a good representative of 
the entire ensemble. 

For completeness let us give here the definition and the relations between thermodynamic observables, i.e. the 
partition sum at temperature T = 1//3, Z(/3,n) = Y,{ n } e ^ f>H({n]) = e~ 0Nf ^\ the Gibbs-Boltzmann distribution 
Mb({«}) = e~^ H ^ n ^ /Z(/3, p), the free energy per particle f(p,(3) = e — pp — Ts, where e = Vdyf, s — —&rf, 
p = —d^f are respectively the average energy and entropy per particle, and the density. 

A. The liquid phase: RS solution 

We now proceed to discuss the cavity method at the replica symmetric (RS) level, extending the derivation of [22] 
to the soft (finite V) model. This calculation provides the exact solution of the problem in the liquid, homogeneous, 
phase. As we mentioned above we consider the thermodynamic limit of random regular graphs with connectivity 
c = k + 1, which locally converge to regular trees. The derivation of the RS cavity equations presented in this section 
starts by considering the solution of the model on a finite tree (a slightly different method will be used in section III B) . 
As represented pictorially in Fig. 1, trees admit naturally a recursive description. Consider the sites 1, . . . , k in absence 
of their common neighbor (they become so-called "cavity sites" ) . Because of the tree structure the Gibbs-Boltzmann 
probability distribution factorizes over the fc-subtrees rooted at 1, . . . , k. Each "cavity site" i can be found in three 
different states: empty if m = 0, saturated if rii = 1 and if the number of occupied neighbors above it (which is k 
at most) is > £ and unsaturated if rii = 1 and if the number of neighboring particles above is < £. We shall denote 
{V>f, 4>f, tpf} the respective probabilities of these three states with respect to the Gibbs-Boltzmann measure. Thanks 
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FIG. 1: Illustration of a cavity iteration for the classical system. The sites 1, • • • , k are described, in absence of site 0, by cavity 
fields ipi' u ' a , I = 1, • • • , k. Integrating out these sites leads to the expression of i/>q'"' s in Eq. (2). 

to the factorization properties of trees one can obtain recursion relations between these cavity probabilities [22, 44]: 
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where Zi tcr is a normalization constant ensuring + V'o + V'o = 1 • Note that the j = term in the sum in the second 
line of the above equation (which is the only one for i — 1) should be interpreted formally as giving a contribution 
equal to 1. Deep inside a very large regular tree, for an homogeneous liquid phase, the cavity probabilities converge 
to the fixed-point solution of Eq. (2), 
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The solution of Eq. (3) can easily be found numerically for any choice of the parameters (3, ^, V. 

The recursive equations written above were strictly valid for trees. The local tree-like character of random graphs 
allows however, under the correlation decay hypothesis of the RS cavity method, to apply these results to the ther- 
modynamic limit of Bcthe lattices. One finds in particular for the free energy density [22, 44]: 



/ — A/ S ite — -A/ii nk 
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- -0 l°g ^sitc + ^ log Z iink 



(4) 



with 



E ( k + x ) m k+1 ~\r + re-p v y + E ( fc + 1 ) (yn k+1 ~\r + re- pv y^ e -^ v 



(5) 



All thermodynamic observables can then be obtained from Eq. (4) by derivation with respect to the parameters 
/?, fi, V. Actually only their explicit dependence has to be derived, the order parameter equations (3) being stationary 
point conditions of / with respect to ip e ' u > 8 . 

The replica symmetric solution provides the natural description for the homogeneous, liquid (or paramagnetic) 
phase where all correlations are short range. However, the correctness of the RS ansatz breaks down when one 
approaches the glass phase, where the system develops amorphous density profiles and long range correlations. In 
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this model this happens at low enough temperature, when the competition between the chemical potential /i and the 
interaction energy V is the strongest. One evidence for the inconsistency of the RS assumption is the fact that the 
associated entropy becomes negative in this region of parameters, which is impossible for a discrete model. Another 
criterion for the incorrectness of the RS ansatz is the divergence of the so-called "spin-glass" susceptibility, defined as 
Xsg(/3, m) = N^ 1 J2i j( n i n j)c- F° r au parameters we investigated we found no divergence of the susceptibility. This 
means that the resolution of the entropy crisis requires a phase transition which has to be discontinuous since it is 
not associated to a diverging susceptibility. This scenario is typical for models displaying 1 step replica symmetry 
breaking (1RSB) [22]. 



B. The glass phase: 1RSB solution 



The replica-symmetric approach stated above was based on an hypothesis of correlation decay, which allowed to 
neglect the presence of long loops in random graphs and to use the results derived on a tree. This assumption fails 
at low temperature, and must then be replaced by the first step of replica symmetry breaking (1RSB). In the latter 
description one partitions the configuration space of the model in pure states, identified by an index a, each of them 
having free-energy f a . One assumes then that the Gibbs-Boltzmann measure restricted to the pure state a enjoys 
the correlation decay property of the RS ansatz, in other words that there is a (non-homogeneous) solution of the 
equations (2) for each pure state a, and that its free-energy f a is given by the RS expression (4). The cavity method 
at the 1RSB level corresponds then to a statistical treatment of these pure states. One introduces the "complexity" 
or "configurational entropy" £(/), such that the number of pure states a with free-energy f a w / is <~ e NT, ^\ The 
computation of this function is done through the partition function [22, 44, 45] 

Z(m;(3,n) = V e _/3JVm/ ° ~ / 4f e JV(E(/) " /3m/) ~ e~^ min ' e tf»in'W](™/-? E (/)) = e -^/W(™;0.AO j (g) 

•//m.n 

where the sum is over all the pure states, and the interval [/ m in,/max] corresponds to the domain where £(/) > 0. 
The parameter m, which is known in the literature as the Parisi parameter, plays the role of a Lagrange parameter 
conjugated to / (an effective temperature). Indeed the "replicated free energy" (p(m) and £ are Legendre transforms 
one of the other [45], 

171 Hm) = mutf^] {mf ~ TE(/)} ' (7) 

and from the properties of the Legendre transform one has 

0/£(/) = mp , S = Pm 2 d m 4>{m) , / = d m [m<f>(m)\ . (8) 

The usual partition function is recovered from Eq. (6) by setting m = 1. However it can happen that the saddle- 
point in (6) occurs for m = 1 outside the authorized interval [/ m in; /max]- In that case the dominant contribution to 
the partition function corresponds to the pure states with free-energy / = / m ; n , which are in sub-exponential number 

as S(/min) = 0. 

The computation of (j>(m) is performed through a statistical study of the solutions of (2) among the various pure 
states. We call rj = {ip e , tp u , -0 s } the set of the (normalized) fields that appear in Eq. (2) and 7^(77) their joint 
probability distribution. Eq. (2) defines a map rjo = 57(771 , • • • , r]k) that allows to construct a new field from a set of k 
fields. Using this map, the 1RSB equations read: 

k 

= \ ff[ dm P(rn) 5{v-v{m,-,m)) z itcr ( m ,..., Vk ) m , (9) 
•* i=i 

where Z is the normalization constant and Zi tcr = exp(— /3A/i ter ), which is defined by the denominator of Eq. (2), 
takes into account the shift in the free energy after a cavity iteration [44]. Once that V{r\) is known it is possible to 
compute all thermodynamic potentials, in particular the "replicated free energy", <p(m; (5 , fi) . In fact, one has: 

c 

<t> = A0 s ite - -A0ii nk (10) 



with 



3m J 



'fc+i 
.1=1 



Z site ( Vl ,...,r, k+1 ) m (11) 



A^i„ k = -J- In / d m d m V{m)V(m) ^iink(?7i^2) m (12) 
pm J 
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and 

/fe+i \ (i-i j ,u + ^ls e -pv 

\i=l / \j=0 l<j 1 <---<i 3 <fe+l p=l Vi p 

fe+i j e -/^\ (13) 

j=£ l<ii<---<i j <fc+l p=l 

Then, from </>(to; /3, /x), using (7) and (8) with to as a parameter, it is possible to compute the curve £(/) for every 
value of (/?, /x) [45]. 

Depending on the value of the external parameters /?, /i, V and m three situations can occur: 

• Only the trivial replica symmetric solution exists at m = 1, which is recovered when V(rj) = 5(t] — r] rs ). In this 
regime the RS approach is correct and the system is in its liquid phase. 

• If at m — 1 there is a non trivial solution of the 1RSB equations with positive complexity the Gibbs measure is 
dominated by an exponential number of states and the system is in its clustered phase. 

• If a non trivial solution exists at m — 1, but the associated complexity is negative, this means that those 
states are absent in the thermodynamic limit. Instead the total free energy is dominated by the states of free 
energy / ma x, that satisfy S(/ m ax) = 0. This condition is equivalent from Eq. (8) to the extremization condition 
d m cf)(m) = 0. We call to* the value of to at which the extremum is attained (note that we are only interested 
in the region to* < 1). In this phase the system is said to be in its condensed phase, as its equilibrium measure 
is supported by a sub-exponential number of pure states. 

The passage between these three regimes defines two transitions which are well-known in the context of glasses (see 
e.g. [46] for a review). The first one represents a dynamical transition at which the phase space of the system gets 
fragmented into an exponential number of states. At this point ergodicity breaking takes place, but the thermody- 
namics of the system is completely blind to the transition [46] . On the other hand, the second transition towards the 
condensed phase corresponds to the truly thermodynamic phase transition. This phase transition, which is usually 
called Kauzmann transition, is second order according to the Ehrenfest classification, but still it has peculiar features 
with respect to usual second order phase transitions, since the order parameter jumps at the transition and the physics 
around it is dominated by nucleation in finite dimension. For these reasons it is also known under the name of Random 
First Order Transition (RFOT) [12]. 



C. Methods of solution of the cavity equations 

The RSB equations (9) can be solved through the population dynamics method [44] , based on the representation of 
the probability distribution "P(ry) as a weighted sample. Starting with an initial representation (population) of V(j]) 
in the form 

M AT 

T^kn) = E Wi s ^ ~ % ) sucil tliat E Wi = 1 ' ( 14 ) 

i=l i=l 

the fixed point of Eq.(9) is found by iterating three main steps: 

1. sample k elements . . . ,rj ik } from V(f]), independently with the probabilities u>i\ 

2. compute a new field r/ = ^(r)^ , ... ,r]i k ) according to Eq.(2); 

3. add to the new population the element r] with weight cx Ziter(%i, ... 7 i] ik ) m . 

The update of the new population can be either parallel or serial. In the first case the same population is used for M 
iterations, after which a new representation is available and it is substituted to the older one. Conversely, in the serial 
update the new generated element is immediately substituted in a random position in the pre-existing population. In 
some cases the parallel update can fail to find a fixed point solution of the self-consistency equation (9) and enters 
instead a periodic orbit (modulation phenomenon) ; this spurious effect can be avoided using the serial update version 
of the population dynamics method. 
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Once convergence is reached all thermodynamic observables can be computed sampling from the distribution. All 
the results of this section have been derived within this scheme where we fixed the size of the population equal to 2 17 
and we checked that finite population size effects were negligible. Both in the classical and in the quantum case we 
used the serial update. The population dynamics technique sketched here is in fact a particular case of the quantum 
population dynamics algorithm, and more details are presented in Sec. Ill C. 

The numerical accuracy of the population dynamics obviously depends on the number AT of representatives of 
the distribution V(i]). However, a large value of N is not sufficient to ensure a good precision of this discretized 
representation: if only a few weights Wi dominate the sum in (14) one should consider that the effective size of the 
population is the number of such dominant elements. To be more quantitative one can define an inverse participation 
ratio 



Ip R=T7 5>H . (15) 




where Wi are the normalized weights of the sampled population defined in (14). Its value is easily evaluated in the 
two limits of a perfectly balanced population and of a single dominant weight: 

and one can thus define an effective population size as Af e g = N x IPR. The inverse population ratio should thus be 
maintained as close as possible to one to achieve a good numerical precision. We observed that the implementation 
of the population dynamics method with the definition of the fields given in Eqs.(2) leads to very small values of the 
IPR (this can be traced back to the presence of the factor e"^ which takes very high values for the relevant values 
of the parameter) . It turns out that a redefinition of the fields allows to bypass this problem. We found much larger 
values of the IPR using 

^iter ^iter '"iter 

where r u and r s are two a priori arbitrary parameters (a good choice for the case c = 3, £ = 1 was r u = 1/5 and 
r s = 2/5). In terms of the transformed fields, the local self consistent equations (given here only for £ = 1 for 
simplicity) become: 



K 

^n(^ + Ce^ r "+^e^) 

iter i=1 



Z' 

iter i=1 



(181 

\i> e< 1 w e n 



Z' 

itor i=l j=l l<ii< — <»^<fc p=l 



itor \i=l »=1 / 

Of course analogous transformations take place in the expression of Z s i te and Zi; n k given in Eq. (13). This change 
of parameters was particularly important in the quantum case discussed below where the population sizes are much 
more limited by computational costs. 



D. Thermodynamics of the classical model 



In Fig. 2 we report the phase diagram for the classical model with c = 3 and I = 1, as a function of density 
and of temperature (in units of V). The two lines denote the dynamic and Kauzmann transitions, below which the 
system is found in the glass phase. Increasing the density from a low value, for a fixed low enough temperature, 
the system encounters first the dynamic transition and then the Kauzmann transition towards the condensed glass 
phase. However if one keeps increasing the density the system returns to a liquid phase; this is because if the chemical 
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FIG. 2: Classical phase diagram in the (p, T) plane for c = 3 and 1 = 1. The dynamical (T<j) and Kauzmann (Tk) temperatures 
are plotted as a function of the density. See Ref. [47] for a study of a realistic model that displays a similar phase diagram. 

potential is large enough (/i > F), almost all sites will be occupied leaving only a small fraction of "holes". The 
system of holes has a short range interaction similar to the original one (i.e. for particles) and is therefore liquid 
when the density of holes is small enough, while it becomes glassy when the density of holes is large enough. The 
"hole glass" transition which is approached from higher density is therefore due to the fact that we introduced a finite 
repulsion V which can be overcome by a large enough p, and it is qualitatively equivalent to the transition which takes 
place at lower densities. A very similar phase diagram has been found in the study of a realistic model for colloidal 
systems [47], which confirms the validity of the Biroli-Mezard model to describe the physics of structural glasses. 

The thermodynamics of the system inside the glass phase has to be computed exploiting the 1RSB formalism 
developed in the previous section. Fig. 3 (a) shows in blue the density as a function of the chemical potential at /3 = 8 
and j3 = 30, both in the liquid and in the glass phase. Then, the red dashed line represents the RS solution within the 
glass region, which clearly differs from the true 1RSB solution. For the same two temperatures the behavior of the 
parameter m* as a function of the density is shown in Fig. 3(b). It differs from one only inside the condensed glass 
phase, after the Kauzmann transition. It decreases fast going deeper into the glass phase, down to a value proportional 
to the temperature. Note that the model presents a glass transition also at zero temperature, as a function of the 
density, with pd — 0.5708 and pk = 0.5725 on the "particle" side [22, 48]. This is of particular interest in view of the 
extension of the model to the quantum case. 

Although we will not be concerned with this physics in the rest of the paper, it is worth mentioning at this point 
that in the limit of zero temperature the parameter m* approaches zero at a finite density Prcp = 0.57574 [22, 48]. 
Below this value of density, the average energy is zero, while above this value it becomes non-zero (the system is so 
dense that it is unable to satisfy the constraint that each particle has less or equal than £ neighbors). This point 
corresponds to what is called "random close packing" or "J-point" in the literature on hard spheres in the continuum 
(see e.g. [26] for a review and a list of relevant references). This is seen in more detail in the inset of Fig. 3(a) which 
reports data for /3 — 30 that are very close to the zero temperature limit: above the glass transition px, the density 
in the 1RSB solution increases slower than in the RS solution and reaches a plateau at the value prcp = 0.57574 for 
fi < V = 1. Only when p > V = 1 the constraints are violated at zero temperature and the density can increase 
above prcp- The latter is marked by the horizontal dashed line in the figure. 

One can recognize that the glass transition of the original model proposed in [21, 22] is recovered here as the 
zero temperature limit of the glass transition, and it is here reached in the "SAT" (namely zero energy) phase. It 
corresponds indeed, to the limit 0V — > oo for a suitable and finite value of ftp. On the other hand, in that context, 
the presence of a "hole" glass transition is excluded a priori by the definition of the model, which does not allow 
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FIG. 3: ie/£ Panel: Average density as a function of the chemical potential, at j3 = 8 (upper panel) and /3 = 30 (lower panel). 
The inset shows a zoom of the /3 = 30 data around the Kauzmann transition. Right Panel: Behavior of the parameter m* as 
a function of the density at /3 = 8 and /3 — 30. 

configurations violating the energy constraint. 



III. THE QUANTUM MODEL 
A. The model 



The quantum version of the model is obtained by adding to its classical potential energy a kinetic term allowing 
hoppings of particles between neighboring sites. Since the classical energy forbids configuration with more than one 
particle per site, we are in fact dealing with interacting hard-core bosons. The Hamiltonian is now an operator acting 
on the Hilbert space spanned by the 2 N classical configurations of the occupation numbers. The potential energy is 
diagonal in this basis. Denoting J the intensity of the quantum tunneling between neighboring sites one obtains the 
following Hamiltonian: 



H = - J ^(ajfij + aja,;) + 



(19) 



•:.i> 



where the first sum runs over the pairs of adjacent sites, a\ and &i are standard hard-core bosonic creation and 

■ t. A thermodynamic study of the model then amounts 



annihilation operators on site i, hi = a\&i and = J2jedi '• 
to the computation of the quantum partition function Z = 



Tr 



and of the average value of the observables O 



defined as (6) = Tr 



-PHQ 



/z. 



As in the classical case we consider the model on the Bethe lattice and we solve it through the quantum cavity 
method, a generalization of the classical cavity method which has been derived first for spin [39, 40], and then for 
bosonic systems [41]. The basic idea of the method is to exploit the Suzuki- Trotter formalism in order to write the 
partition function as a path integral. The Hamiltonian is naturally split into two non-commuting term H = Hi + H2, 
with Hi the kinetic energy and H2 the classical potential energy. The Suzuki- Trotter method amounts to work in the 
diagonal basis of one of the two terms of the Hamiltonian, and treat the non-commutativity of the two terms through 
the Suzuki- Trotter formula 



Z = Tr 



-f)H 



— lim 

]V„->-oc 



Tr 



-H-2 



(20) 
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The computation which leads to the path integral formulation proceeds with the insertion of N s copies of the identity 
between the N s factors of this product. In particular for the purpose of the implementation of the cavity method, it 
turns out that the most convenient basis to use in this operation is that of the occupation numbers. This is due to 
the discreteness of this basis. Then the quantum system is mapped into a classical one with an additional "imaginary 
time dimension". This procedure, which preserves the spatial structure of the graph of the interactions, allows to 
apply the general formalism of the classical cavity method. The crucial difference is that one must now consider as 
basic variables the whole trajectories in time [0, f3] which describe the state on each site. Detailed explanations of the 
bosonic quantum cavity method can be found in [41] which studied the Bosc-Hubbard model on the Bethe lattice. In 
the following we shall therefore focus on the new features that arise here due to the specific form of the interaction 
(1) between neighboring sites, which was absent from the treatment in [41] (see however Ref. [37] for a study of the 
extended Hubbard model, that corresponds to the present model in the special case £ = 0, but has a very different 
behavior from the I > in which we are interested here). The next subsection (Sec. IIIB) will present a derivation 
of the quantum cavity equations which describe this system, while the numerical method to solve them will be the 
subject of sections IIIC and HID. This part of the paper is rather technical and the reader who is not interested in 
the details of the derivation can safely jump to Sec. IV where we shall discuss the results of this procedure and the 
physical properties of the model. 



B. Derivation of the cavity equations 

1. Rewriting the partition function as a simple factor graph 

We use in this section the same notation as in [41] for the imaginary time copies of the occupation variables. We 
consider a discrete number of Suzuki- Trotter slices labeled by a = 1, • • • , N s . We denote by nf the occupancy of site 
i at Suzuki- Trotter time a, by rii — (nj, • • • , nf s ) the imaginary time trajectory of site i, by n" = (n", • • • , n^) the 
configuration at time a, by n = (n 1 , ■ ■ ■ ,n N ") the imaginary time trajectory of the full system. 

We write the partition function 



Tr 



— lim 

Af B -s-oo 



E ex p 



^EE^w-K^}) 



Ns N 



K 



a— 1 i—l 



a=l 



J2 J(a\&j+a^ai) 



1 



n a+ ) (21) 



where 



«(»?. Ke*}) =Vnf q? %?) - fjmf 

9? = E n i - 1 



(22) 



and we write here and in the following 9(x) = l(x > 0) and 1(A) = 1 if the condition A is fulfilled, otherwise. It is 
simple to show that [41] 



(jf\e 



Yl J(a\aj+ajai) 



i >=sn 



ni(< +i =<+E[^-^ 

j£di 



— s-J J 



o 



(23) 



where we introduced variables y with the following notations: we denote by yf_+j € {0, 1} the variable that indicates 
an hopping event from site i to j at time a; by = {yfUj,yf^i} the two hopping variables on link (i, j); and as for 
occupation variables, we use bold notation to indicate the imaginary time trajectory of a variable and an underline to 
indicate the full configuration of the y's. Finally, we take by convention x y=a — 1 for any value of x (including zero). 

Eq. (23) can be checked by inspecting the behavior of its left and right-hand-side order by order in 1/N S . The 
leading term corresponds to n a = and indeed all y's must vanish at this order. We note that the partition 

function has the following form: 



N 



Z = J [m E II w h(y(t,j))T[ w s(ni,{ n j}jedi,{y{i,j)}jedi) 

iv a — >oo 



(24) 
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with 



w h 



(v<ij>) = n 



a=l 



N s 



(25) 



and 



w s (ni,{nj}j e di,{y(i,j)}jedi) = exp 







exp 


N s 






n< 

Q = l 


!■(• 






exp 





]>>«,K e9 J) 



^s(«*){j/(i,j)b€8i) • 



(26) 



One can easily recognize that the graph of the interactions in the above representation has small loops, so, at first 
glance, it may appear not suited for cavity calculations. In order to avoid these loops and thus to set down an 
appropriate cavity treatment of the problem we adopt a trick which consists in copying on each site the variables of 
the neighboring sites and the jumps across the edges connected to the site. Let us introduce then 



^ = {nf,{nf} jeBi } 



(27) 



where {fi} a }je8» are the copies on site i of the occupation numbers of its neighbors at time a. In addition, we copy 

each variable yf^j on the two sites introducing Uj%i = v]%i = Uj-t-n an d on eacn site we store these variables in 
a variable T" : 

T = {y\Z j} hedi = {y&j, y&ibew • (28) 
Then we impose the consistency among neighboring sites' jumps and occupation numbers through a constraint 



Wi 



(is;, Y,-; i/.-, T.) = TT S a jc,5„ a „i a S ia ja 5 ia ,« 



(29) 



a=l 



Putting together this constraint and the hopping weight Wh{yu,j\) (properly expressed in terms of the new variables 
T) we obtain a total weight on each link 



wiiv^T^Vj.Tj) = w h {'T i ,T :j )wi{v i ,T i ;i'j,Tj 
The partition function now becomes 

N 

^ TT wdui, T, : ; i/,J,) 

Ar B -s-oo 



Z= lim y^T\wi(vi,'Ti-,Vj,'Tj)T\w s {vi,~£i). 



(30) 



(31) 



i=l 



where we have defined Y = {Yj} and the argument of w s has been expressed in terms of the copied variables {i/j, Yj}, 
living on site i. 

The last form of Z is a standard form with local weights and link weights only, which is therefore defined on the 
original graph and hence without short loops; it is therefore suitable for a cavity treatment and the cavity equations 
as well as the Bethe free energy can be written straightforwardly. Introducing fji^j (i/j, Yj) = Z^j(yi, Y^/z^j, 
which represents the marginal probability of the variables fj and Yj defined on site i, in absence of site j, the cavity 
equations read [44] 



f]i^j(vi, Yj) = — — w s {vi, Yj) ] [ y~] fjk^ijvk, Y k )wi(uk, Tfc;i/j, Yj) . 



(32) 



kedi\j v k ,T k 



12 



The Bcthc free energy associated to the partition function, expressed in terms of the copied variables as in Eq.(31), is 
1 1 N ( \ 

1 1 i=l \v % ,T % j Hi v T J 



(33) 



This expression is exact whenever the underlying graph is a tree while it corresponds to the Bethc approximation for 
a general graph. 

2. Simplification of the cavity equation 

The cavity equations written above are not very practical to handle, since there is a lot of redundancy in the copied 
variables. We can eliminate much of this redundancy by using the delta functions in the weights iwj. 
First we make explicit the dependence of the cavity field on the original variables in the following form: 

^•(i/j.Ti) = jh-fj(n i ,n5,{n'fc} fceaiVi ,y| iiJ - ) ,{y| iifc) } fce 8 iNj -) . (34) 
Next we introduce a new cavity field: 

Vk^i{nk,n^,{nf} ledk \i,y^ k) ) = w h (y k {ik) ) ^ %-»(™fc, n k , {nf }z G(9fc \ 4 , y\ %k)l {y^ k) }i e dk\i) , (35) 

and we note that using the delta functions 

E %^i(i/k,T fe )u;j(i/ fc , Y^i/*, Y*) = ^ Vk^i(n k , n t , {nf }i e9 k\i, y\i,k)) ■ ( 36 ) 

This allows to write closed equations for the new cavity fields as follows: 

»?i-^j(*»i."j-.{»»fc}fc€8i\j'V<iJ>) = 7^— J>) 

^2 w s (n»,n5JnUfcea»\j)l/(i,j>){l/(i,fc>}fe€8i\i) II E %^i( n fe' n i'{ n f}ieafc\*'f(i,fc>) ■ ( 37 ) 

{l»<<,*>}keei\j feea»\j {n*}, ee *\< 

We notice at this point that the upper indices of the variables in the equation above are redundant: they can all be 
eliminated by renaming the variables. We obtain the following equation: 

Vi^j(ni,nj,{n k } ke di\j,y{i,j)) = -^—w h (y{ij)) 

^ Ws{ni,nj,{n k }kedi\j,y(,i,j),{y(,i,k)}kedi\j) JJ E n *> { n j}«eafc\*)l/(»,fc>) ■ ( 38 ) 

{W<i,k) }fceai\j kedi\j {ni}i eakV 

To further simplify the equations, we introduce a variable s?^ • which indicates if a site i at time a is saturated by all 
its neighbors but j 

€+i = E < - e > S U = > (39) 

kedi\j 

where we recall that we assume 0(0) — 1, then we consider 

T] i -, j (n i ,n j ,s i -, j ,y {iJ) ) = ^ Vi->-j(ni,nj ,{ n k}kedi\j,y(i,j)) Sa^fi^^) (40) 

{n-k}k£Oi\j 
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with 5 s ^ ]fi{q .^.) = J! 5 sf^,e(gf^) ■ Using the relation 

a— 1 



which can be checked by inspection of the cases q^j ^ — 1, Eq.(38) becomes 



(41) 



f]i^ j (ni,n j ,s i ^ j ,y { ij)) = -^—w h (y {iJ) ) ^ w s (rii,y {itj) , {y(i,k)}kedi\j) 

\V{i,k) \kedi\j 

J2 e-^^<( y ^K+^)-^)^^ e(q ^ ) 11 Vk^(n k ,n u s k ^,y {hk) ) . (42) 

{nk}keai\j kedi\j {s k ^i} 

Now we note that the function r]i- > .j(ni,nj,Si- > .j,y^j^) has a dependence on the variable n 3 - of the 

__/3_ v^^s v a a a 

form r)i^j(ni,nj,Si^j,y^j)) = e jv b 2^=1 n < s i^ n j fji^j(ni,Si^.j,y^j)), so, in terms of the marginal 
»?»->j(n*,Si->j,l/(i,j>) we finally obtain 

- /■ \ w h(y{ij)) \ - f > n . , x 

r)i^j{ni,Si-n,y(ij)) = — - — ■ — 2^ r WiteT{n i ,s i - l .j,y^ j ),{nk,Sk^i,y{i,k)}kedi\j) j [ rj k ^i{n kl s k ^i,y {i ^) 



c t->3 



{ n k}k€di\j 
{V(i,k) }k£di\j 



k£di\j 



where 



Witer(n», Si->j, y(i,j), { n k, Sk->i,U{i,k)}kedi\j) = w s (rii, y{ij),{y{i,k)}kedi\j) 

k£di\j 



e -WJ Sail ^™?™fc s fc^i 



(43) 



(44) 



From Eq.(43) we see that s£L^ enters in the equation only through its product with nf. This means that in order to 
characterize the state of each site one has actually to consider the variables {nf , nf sf^j} and sum over sf^j in case the 
site is empty. Roughly speaking, defining e = {0, 0}, u = {1, 0} and s = {1, 1} we can interpret fji^ 3 {rii, Sj_>.j, yu,j\) 
as the marginal probability over trajectories of the kind [0,/3] — > {e,w, s} which describe the state of each "cavity 
site". Finally one can note that the variable s^j is completely determined by the neighbors' occupation numbers. 
This ensures that, even if locally in order to write the cavity recursions one needs three possible states to describe 
each site, globally the Hilbert space has size 2 N . 

In terms of the cavity marginal probabilities fji^j^m, s^j, y{ij)) the occupation site trajectory on site i is then 
expressed as 



(45) 



{y<.,j)}jea. 



which allows to compute all the local observables. From (33) we can recover the free energy as a function of the last 
defined messages: 



N 



IE* 



\ 



E Mm,{nj},{V(i,j)}) II S^->*( n i' a J-^'W<M)) e ~. E -= inf 

ni ,{nj} jeai jedi s i^i 
\{V(i,j) }jedi 

( \ 



E 



V(i,o) 



e 



~Vi—>-j iP'ii s i^jj V{i,j) )%— H ( n j ) s j— >i j U{i,j)) 



(46) 



/ 
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In the case of a regular Bethe lattice where all sites have degree c, denoting ?7 cav the common value of the distributions 
fli-tj on all edges, we obtain for the RS free energy density: 



* = -m [nZ 



( 

w s (n ,{nj,y j}) jf| ^ Vca,v{nj, Sj^ , Voj)e~ 

«o je{i : ..., c } s i^o 
\{ n j>yoj}j<={i,...,o} 

( 



\ 



3V v^N s a a a 



+ — In 

2/3 



E 



e N * 



777 Soil "o n i ( s o^i+ s i-+o) 



I ™l,*l->0 



wh(y) 



(n , s -n, y)%av(«i) sno, y) 



(47) 



C. Resolution of the cavity equation for the random regular graph 

In this section we explain the idea and the method used to solve the quantum cavity equation for a random regular 
graph. This task is rather demanding since already at the RS level it consists in a functional self-consistent equation. 
However, as it was shown in [40, 41], there is a representation of the equation which allows for its numerical resolution. 

Let us define hi = {rii, V{ij)} the whole set of variables argument of the cavity marginals in Eq.(43), which 
determine the "field" on the i site when j is absent. In the following we deal with a random regular graph, so for 
simplicity, we will label with the root site which is added in the cavity iteration when we merge k = c — 1 neighboring 
sites and their associated fields {^i}ie{i,...,fe}- Finally, we label with an index k + 1 the new site which is missing 
"downwards" with respect to the root. The calculation is based on the observation [40, 41] that for a random regular 
graph Eq. (43) has the form 

rj° cav (h ) = — J2 Z(h u ...,h k )P(h \h 1 ---h k )rj 1 cav (h 1 )--- V k cav (h k ) , (48) 
hi—hk 

where P is a conditional probability, positive for all arguments and normalized and 

Z(hi, . . . ,h k ) = ^2 Wh(y(o,k+i))witer{no,s ^k+i,y(o,k+i),{ h j}je{i,...,k}) ■ ( 49 ) 

/io={no,So_).fc+i,!/<o,js+i} } 

We note that at the RS level we look for the homogeneous solution so all r^ av = rj cav are equivalent, however in 
equation (48) we keep the general form to allow the generalization to the 1RSB case. 

As detailed in [40, 41], a strategy of resolution is based on the representation of 77(h) as a weighted sample of 
trajectories: 

AA traj 

?7cav(h) = E 9i S ( h - h i) ) ( 50 ) 

i=l 

where the A/t ra j weights of the trajectories are normalized according to 

E 9i = 1 ■ (51) 

One should not confuse this population of quantum trajectories (that represents the RS solution of the quantum 
problem) with the population of cavity fields that we introduced in the classical case in Eq. (14) to describe the 1RSB 
solution. The form (50) provides an approximate representation of r] cav which becomes more and more accurate as 
■A/traj grows. Starting from an initial weighted sample of 77 cav in the form of (50) the procedure amounts to iterate the 
following three steps: 

1. extract k trajectories hi = {rij, s^j , yuj) } in [1, . . . , A/traj] from (50) according to their weight gi 

2. draw the new trajectory from P(ho\h\ ■ ■ ■ h k ) 

3. set h\ — h and g[ = Z(h\, . . . , h k ) 
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where the second point amounts to a single site problem that will be further explained in the next section. Once a 
new representation of ry cav in terms of h[ is available it can be substituted in the r.h.s. of Eq.(48) and the procedure 
is iterated until the convergence, in the sense of the observable expectation values, is reached. 

Before entering into the description of the generation of the new trajectory given the k neighbors, we anticipate 
that, as in the classical case, also the quantum model is characterized by a glass phase at high enough densities, where 
the RS equations are no more correct (they have the same pathologies as in the classical case, yielding for instance 
negative entropies). The 1RSB treatment is needed also for the quantum problem and within this framework we have 
derived all the results concerning the glass phase, presented in the following. 

At this level however, the generalization of the classical 1RSB equation (9) is straightforward. The only difference 
is that now the fields h — {ip e , ip u , ip s } of Eq.(9) are probability distributions over paths, namely the "cavity fields" 
?7cav 01 the RS equations (48) and Z- lter is the normalization constant z cav that appears in the same Eq.(48). Then, 
the 1-RSB cavity equation for the random regular graph reads: 

Pfacav] oc J dP^] . . . dP[ v k cav }S ( T]c&v fijlcav) ' • • i ^7cav)) ^cav i 

(52) 

where f{r]l av , . . . , ?7cav) is the r.h.s. of Eq.(48), and a normalization constant has been hidden in the notation oc. The 
numerical resolution of the equation (52) necessarily goes through the representation of P[?7 cav ] as a population of 
"fields" 7? cav 

N 

■P[?7cav] = X] W i S (Vc a v - Vi) . ( 53 ) 
i=l 

where the weights derive from the cavity iteration w oc z™ v (analogously as in the classical case) and each r] cav is itself 
a population of trajectories. This procedure is of course accompanied by a demanding computational cost which is 
however still attainable and controllable. For the interested reader we refer to [49] for more details about the 1RSB 
cavity method in quantum systems and the computation of observables. 

D. Generation of a trajectory 

In this section we will focus on the problem of the generation of trajectories for the thermodynamic limit of a 
random regular graph of connectivity c. As we said before the strategy of resolution is based on the representation 
of rj(h) as a weighted sample of trajectories. Then one will pick k trajectories hi ■ ■ ■ hk from the corresponding 
distributions, and use P{ho\h\ ■••hk) to draw a new trajectory ho, to which a weight Z(h\, . . . ,hk) is assigned. 

In order to calculate Z and to generate a new path we need to take into account all events encoded in the neighboring 
trajectories. In particular we can distinguish between three of them: 

• a particle jump from or towards the neighboring sites, which is encoded in k} 

• a change in the {s»->o}»e{i,...,k} variables 

• a jump from or towards the neighboring sites along the edge connecting the site under consideration, described 
by {y<i,o>}ie{i,...,fc} variables. 

The last class of events is a subset of the first one, but it is important to take it into account separately because 
it induces constraints on the new occupation number trajectory. All the other events, instead, will enter as shifts 
of the effective chemical potential induced on the added site. Moreover we stress that in the continuous time limit 
(N a — > oo ) the hopping trajectories h typically contain only a finite number (with respect to N s ) of events. We can 
thus assume that they occur at different values of the discrete time for the different trajectories. We call p the total 
number of hopping events occuring in (hi, . . . , hk)- We also denote a\ < ■ ■ ■ < a p their discrete time of occurence. 

We consider the Hilbert space of a single site, with a and a) the annihilation/creation operators. We introduce the 
operator bj — a (resp. bj = at) if at time aj there is a jump towards (resp. outside) the site under consideration, 
while bj = 1 when there is an event which does not involve the new vertex, but only the neighboring trajectories. 
The former case corresponds to the occurrence of an event of the third kind, while the latter corresponds to an event 
of the first or second kind. Finally we introduce c a — bj when a — aj, c a — 1 (the identity operator) otherwise. In 
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terms of these operators we get the following expression for Z: 

Z(h X ,. ■ ■ ,/lfe) = Wh(h )Witei(h , {/»»}»£{!,.. 



ho 



(54) 

= 2 II «|e^ (Aia({<,s ^ o}ie{1 - - fc})Qta+J(a+af)) c Q K +1 ) 

no a=l 

up to corrections of order N s ~ 2 , where we have introduced an effective chemical potential, which depends on the state 
of the neighboring sites and is defined as follows: 



M a (K,*f\ }i€e{i,...,fc}) = -V (c+ fe+ i%o-+ fe+ i) +J2< s ^oj +/" 



(55) 



To take the continuous time limit it is convenient to define tj — -fi-ctj, which are the continuous times of the 
hopping events in (hi, . . . , hk), and to denote /Uj the common value of the fi a for a e [ctj-i, a,]. We also introduce 
Wi(X) = e A (Mia t a+^(a+a t )) ; the propagator of an imaginary time evolution on an interval of length A for a single site 
Hamiltonian Hi = —[liC^a — J(a + a 1 *). This propagator is a two by two matrix and it easily allows to compute Z in 
the continuous limit according to the relation: 

Z(h u . . . ,h k ) = Tr (w^rO&iWbfa - ti)& 2 ...W p (t p - Tp - 1 )b p W p+1 (fi - T P )) . (56) 

We can now look at the process of the generation of a new trajectory ho given the ones of the c — 1 other neighbors, 
which respects the following probability law 

, > w h (h )w ite Aho, {fci}ie{i,...,fc}) ,.-s 

P(h \h 1 ,...,hk) = ^tt 7—r 1 — ■ (57) 

Z(hi, hk) 

The general scheme to determine h = { n o, s o, y<o,fc+i) } consists in first drawing the occupation number trajectory 
no and then deduce y/o,k+i) from its jumps not associated to the events in (y/o,i) 5 • • • , y<o,k))- Finally the trajectory 
of s is completely determined by the neighboring ones (n\, . .. ,rik) and does not have to be generated. We keep the 
notation t\ <■■■ <t p and b\, . . . , b p for the continuous time of the events in (hi, . . . , hk)- Let us call n = n(r = 0), 
rii (resp. n'j) the value of n(r) at a time just after Tj (resp. just before Tj+i), with the conventions n' p = no- The joint 
probability law of these occupation numbers which arises from the expressions of Wh and u>iter given in Eqs. (25) and 
(44) reads in the continuous time limit 

1 p 

P(n ,n ,m,ni, . . .,n p |hi, ...,h k ) = -— j^r{no\Wi(Ti)\n ) J[ {(n'^b^ni) (ni\W i+ i(T i+ i - r 4 )|^)| , (58) 



with Tp+i = (3. This probability law is well normalized according to the above expression of Z(hi, . . . , hk). It follows 
immediately from the above equation that whenever the operator bi is non trivial, i.e. when there is a jump along 
the edge connecting the neighbors to the new site, the state of the site is uniquely defined by the direction of the 
jump. This is a trivial consequence of the fact that we are dealing with hard-core bosons. Then, we have to generate 
only those values in the sequence (no, n' Q , m, n[, . . . , n p ) which are not already fixed. This is a task that can be done 
sequentially, starting from no, quite easily. Once all the intermediate occupation numbers are fixed, the generation of 
the rest of the path will proceed independently for every interval {r^ + i — Tj}. Inside each of them, in fact, one has 
to generate a trajectory with fixed boundary conditions according to the path integral representation of an effective 
Hamiltonian which depends on time interval Hi = —Hiata — J(a + a^). All the details for this procedure are explained 
in [41] for the Bose-Hubbard model. At this level, in fact, the task is completely equivalent to that of the Bosc- 
Hubbard, despite the interaction between neighboring sites, which is here hidden in effective single site Hamiltonians. 
In the next section we will give more details on the parameters we used for the computations presented in this paper. 



IV. RESULTS FOR THE QUANTUM MODEL 



The system investigated here presents a rich phase diagram, which emerges from the classical case discussed in 
section II when the additional effect of the hopping is taken into account. In this section we first discuss the finite 
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temperature phase diagram, which is directly accessible to our method, and then we argue on how it extrapolates to the 
zero temperature limit. We also present exact diagonalization data for the ground state to support the extrapolation 
of the cavity data to zero temperature. 

Note that the hopping amplitude J plays in this model the role of the quantum parameter T = A/a (the ratio of 
the De Broglic thermal wavelength and the interparticle distance) defined in the introduction. As we showed in the 
previous sections, for the classical problem where J = 0, the thermodynamic glass phase is delimited by the curve 
Tk{p) of Fig. 2 (while the dynamical glass is delimited by Td{p)) in the temperature-density plane, similarly to what 
happens for more realistic models of structural glasses [47] . One can then imagine to approach the glass phase from 
higher or lower densities. As we discussed in section II D, the two transitions are qualitatively equivalent and in the 
following we will focus on the region of the phase diagram corresponding to small densities. From a physical point of 
view, this side of the glass transition is the most interesting one, since it corresponds to a packing problem, which we 
have relaxed through the soft constraint V. From now on we will focus on the model defined by c = 3 and i = 1. As 
before we will measure energies in units of V throughout this section, leaving it implicit in the text (but not in the 
figures) . 

A. Order parameters 

1. Edwards- Anderson parameter 

As we already discussed for the classical case, a fingerprint of the glass transition is the appearance of a local 
inhomogeneous density profile, signaling the breaking of translation invariance. If one considers the Gibbs measure 
over imaginary-time paths which is constructed by the Suzuki- Trotter decomposition, the transition that separates 
the liquid phase from the glass is of the same kind of the one described in section II B for the classical system. 
It happens in two steps: at a first transition {dynamical transition), the Gibbs measure becomes dominated by an 
exponential number of states, which are counted by a well defined thermodynamic potential, the complexity. As we 
said the thermodynamics of the system is not affected by this transition, while the dynamics is. This has been well 
established in the classical case [46]. In the quantum case, it has been studied in [29, 30] for the imaginary time 
dynamics and in [50, 51] for the real time dynamics. An heuristic argument that relates the transition in the real and 
imaginary time dynamics has been proposed in [9]. At a second transition (Kauzmann transition), the complexity 
vanishes, the system undergoes the true thermodynamic glass transition to the glass phase which is characterized by 
a sub-exponential (in the size of the system) number of dominant states. 

As in section II B, we can label the different glass states by an index a; each state has free energy f a and we can 
perform the same construction as in section II B to compute the complexity for the Gibbs measure on imaginary- 
time paths. We denote by (0) a the average of an observable O in the restriction of the Gibbs measure to a given 
glass state. Inside each glass state one has (hi) a ^ jjJ2i( n i) — Pi while the average density is the same for each 
state, p = jfJ2i ( n i) a - These density fluctuations are conveniently characterized by the Edwards- Anderson order 
parameter [36]: 

^^atEE^K^-p] 2 , (59) 

i a 

where we take the average over the states according to their weights W a cx cxp(— f3N f a ), and the over bar indicates 
the average over random graphs. The order parameter qEA jumps to a finite value at the dynamical transition, where 
the glassy states appear for the first time [46]. This non-zero value of qEA corresponds physically to the long time 
limit of the density-density correlation function [46], as we will discuss later on. As we already said, the jump of the 
order parameter is a peculiarity of the glass transition which still remains second order from the point of view of the 
singularity of the free energy. Note that one should not confuse qEA with the following quantity: 

<Zav = ^ -PF^EE ^7 (fa>a - P)«^) 7 ~ P) > ( 60 ) 

i i a, 7 

The latter quantity is non-zero at finite N, but it should vanish in the thermodynamic limit, since in that limit the 
full Gibbs measure (i.e. the average over all glass states with their weights) must become uniform and (hi) — > p for 
all i. 
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2. Condensate density 

In presence of a finite hopping strength, the system might display Bose- Einstein condensation (BEC). This is 
manifested by a finite expectation value of the bosonic operator (a) ^ in the grand-canonical ensemble. However, 
since the number of particles is conserved by the Hamiltonian and by the cavity equation, in absence of symmetry 
breaking terms the expectation value (a) is always zero. More precisely, one possibility to define the order parameter 
of the BEC transition amounts to add a small perturbation to the Hamiltonian which breaks that symmetry Hh = 
H — hJ2i(<ii + &i) an d then send it to zero after the thermodynamic limit is taken 

^ = l [Ul n J™ Wfc . ( 61 ) 

where (•) h denotes the thermodynamic average in presence of the perturbation. Note that this procedure does not 
present any particular difficulty for the cavity method, that gives direct access to the thermodynamic limit. However, 
an equivalent and more practical way to allow the symmetry breaking, that we used here, is to initialize the population 
dynamics algorithm with a fraction of "asymmetric trajectories" , i.e. trajectories such that the number of hoppings 
from a given site i to a given neighbor j is not equal to the hoppings from j to i. In absence of BEC, these trajectories 
disappear under the cavity iterations, while in the BEC phase the symmetry breaking is preserved, thus allowing 
for a non zero value of (a). The mathematical reasons for this have been detailed in [41]. The physical intuition is 
the following: the cavity field represents the effect on a given site, of the part of the graph that has already been 
integrated out. In the insulating phase, the hoppings are localized, therefore loops are irrelevant and if a boson hops 
off a given site towards one subtree connected to it, it must at some point come back from the same subtree. On the 
other hand, in the BEC phase, there are long range exchanges, therefore a boson can hop off a site to a subtree, travel 
around the (infinite) loops, and come back from another subtree. The "incoming" and "outgoing" currents from a 
site to a given neighbor do not balance (of course, they balance on average, but not on each individual trajectory of 
the population). Indeed, defining a notion of "superfluidity" on the Bethe lattice is tricky, and the idea of bosons 
hopping through infinite loops is probably the closest one to the notion of "winding numbers" on a regular finite 
dimensional lattice [52] . In the following we will refer to the BEC phase as "superfluid" , but one should keep in mind 
this difficulty. 

The order parameter (a) defined above also signals the presence of off diagonal long range order in the system [53] . 
In a uniform phase such as the liquid phase, one can define the condensate density as [53] 

p c = lim (al^ = WF- (62) 

If the latter is finite, there is off diagonal long range order; p c represents the number of bosons that are condensed in 
the ground state divided by the volume. In a non-uniform phase, one has 

i 

for each glass state, and one can define the average over states and graphs: 

P C = ]TW^. (64) 

a 

In a non-uniform phase, the relation with the one-body density matrix is more complicated and we will discuss it later 
on. For the moment, we do not need the relation with the one-body density matrix, since the cavity method gives 
direct access to the expectation of (di) a , and then to the condensate density defined in (64), through the procedure 
described above. 



B. Details on the numerical resolution of the cavity equations: 
results for the order parameters, and how to detect the phase transitions 

In this section we give the details of our implementation of the procedure described in section III to solve the cavity 
equations and investigate the properties of the different phases. Our code is available upon request. In the following, 
we will give indicative values for the different parameters such as population sizes, number of cavity iterations needed 
to equilibrate the populations, etc. These parameters are not fixed a priori but changed slightly depending on the 
region of the phase diagram under investigation. As far as possible given our computational resources, the results have 
been checked to be robust with respect to these parameters by increasing the population sizes in some representative 
points. 
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FIG. 4: Left Panel: Order parameters of the superfluid (p c /p, black circles, left scale, RS cavity method) and glass (c[ea, red 
squares, right scale, 1RSB cavity method) phases for c = 3 and 1 = 1, for three representative values of the hopping J at the 
same temperature j3V = 15, as functions of the chemical potential p. The first order transitions are accompanied by hysteresis; 
full lines corresponds to following the solution upon increasing p, dashed lines correspond to decreasing p. Arrows mark the 
first order transition points p c as determined by the Maxwell construction. See text for a more detailed description. Right 
Panel: Sketches of the free energy as a function of p for the corresponding left panels: from left to right, black lines represent 
the superfluid, red lines the normal liquid, blue lines the glass; p c and pk indicate the phase transitions (we do not report data 
because the jump in the derivative of / at the first order transition is very small and would be invisible in the figure). 



1. Replica symmetric calculations: the superfluid-normal liquid transition 



We started by performing replica symmetric calculations using a population of ~ 32000 trajectories. The population 
was initialized at some low density, with most of the trajectories having no jumps, and a small fraction of trajectories 
having two jumps towards different sites in order to break the symmetry and allow a superfluid phase, as discussed 
above. We then performed ~ 100 cavity iterations to find the fixed point of the population, and then we performed 
~ 500 additional iteration to collect statistics on the averages of the observables. We fixed J and T, and we increased 
p in small steps, in order to reach higher densities, repeating the procedure above but starting from the converged 
population of the previous step (we inject a small fraction of trajectories that break the symmetry, again in order 
to allow a superfluid phase). In Fig. 4 we show the behavior of the order parameter of the superfluid phase, p c f 1 p 
(black circles), as a function of increasing p (corresponding to increasing density). We observe two different behaviors 
depending on the hopping strength: 

• At high temperature or small hopping (small j3J), we observe that p c /p vanishes continuously at some p c 
corresponding to the superfluid transition (upper panel in Fig. 4), which is then a second order transition. The 
density is also continuous at the transition. Close to p c , we observe a power-law behavior (a) ~ \p — pc] 1 ^ 2 , 
with an exponent 1/2 typical of the mean field nature of the underlying lattice; correspondingly, the condensate 
fraction vanishes linearly. 

• At low temperature or large hopping (large (3 J) we observe that p c j p jumps abruptly to zero at some value p c , 
indicating that the transition is first order (middle panel in Fig. 4). In this case, we observe hysteresis when 
we perform increasing or decreasing p scans (see the middle panel of Fig. 4). As a consequence of the first 
order phase transition in the p variables, the density, which is defined as p — —df(p)/dp, has a jump at p c , 
which implies the existence of a region of phase coexistence for p G [ J o_, j o + ] in the canonical ensemble. The 
determination of p c and of p_, p+ can be done in two equivalent ways: either by looking at the grand-canonical 
free energy as a function of p (which is directly accessible within the cavity method) and determining the point 
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at which the free energies of the two phases cross, or Legendre transforming the free energy as a canonical 
function of p and then performing the Maxwell construction over it. Since the jump in density is extremely 
small at the transition, p + — p_ <C (p+ + p_)/2, we found that the Maxwell construction is numerically more 
precise than looking directly at the slope of the free energy —df(p)/dp, which is dominated by the average 
density. 

2. 1RSB calculations: the glass transition 

Next, we turn to the 1RSB equations to investigate the stability of these phases with respect to glassy order. The 
1RSB equations are formulated in terms of a population of populations of trajectory. We typically used 1000 popu- 
lations, each made of 6000 trajectories, which is already extremely demanding in terms of computational resources. 
Fig. 4 shows the behavior of the order parameter of the glass phase, qea (red squares), obtained through the 1RSB 
cavity method, as a function of the chemical potential. In this case we initialize the population at high density (where 
we expect a glass phase) in the following way: almost all trajectories are classical (no jumps), with a small fraction 
breaking the symmetry as usual, and we initialize each population of trajectories at a different average density in 
order to start with a finite qEA- We then perform ~ 150 cavity iterations to find the fixed point and <~ 700 iterations 
to collect data. Since in this case the code is very slow (it can take up to 3 days to compute one state point), in 
order to speed up the computations we run the code independently for each J, T, p with the initialization described 
above. Moreover for each choice of J, T, p wc run ~ 5 computations at different values of m in order to find the 
maximum m* of the free energy. In some cases we performed scans by decreasing p starting at each time by the 
previously converged population: we did this in order to observe hysteresis (dashed red line and squares in the lower 
panel of Fig. 4). Generically, we observe as expected that m* increases on decreasing p, until it reaches m* = 1 at 
the Kauzmann transition px- For lower p, one has m* — 1 and by further decreasing p the non-trivial 1RSB solution 
is lost at the dynamical transition pd, where qea jumps to zero and one gets back the RS solution. In general, we 
found that pk and pd are always extremely close to each other, as in the classical case. A precise determination of 
Pk can be done by looking at the point where the complexity at m = 1 vanishes. On the contrary, the determination 
of pd is much more difficult since it corresponds to the point where the non-trivial 1RSB solution disappears and 
qEA = 0; this is a kind of "spinodal point" and its precise location is very sensible to details of the computation such 
as population size, initialization, number of iterations, etc. Although we managed to get reliable data for pd, in the 
following we focus only on pk, keeping in mind that the two are very close and follow similar trends. Finally, it is 
worth to note at this point that we did not find any superfluid non-trivial 1RSB solution: we comment on this in 
more details in the following. 

3. Summary 

To summarize, we observed generically three different behaviors upon varying p, that are summarized in the three 
panels of Fig. 4 (we set (3 = 15 in the figure): 

1. The upper panel, at small hopping (J = 0.15), represents a second order superfluid-normal liquid phase transi- 
tion, followed by a liquid-glass transition at higher p. The system is superfluid until p c vanishes at p c , then it 
is a normal liquid in an interval of p, until the parameter qEA jumps from zero to a finite value signaling the 
dynamic glass transition pd- Upon further increasing p, the complexity vanishes at px and the system enters 
the glass phase, where m* < 1 (note that pd and pk are indistinguishable on the scale of the figure). Both 
the transitions are second order in the sense of Ehrcnfcst classification: the first derivative of the free energy is 
always continuous and the second derivative is discontinuous at both transition points. 

2. The middle panel (J = 0.2) shows a first order transition between the superfluid and a normal fluid, followed 
again by a liquid-glass transition. Here the condensate fraction p c / p jumps suddenly to zero at p c . In this 
regime the free energy of the superfluid intersects that of the normal liquid at p c with a discontinuity in the 
first derivative. As in the previous case, the glass transition happens here at a higher value of p where qEA 
jumps to a finite value and the free energy of the glass grows smoothly from that of the liquid with continuous 
first derivative. The figure also shows the hysteresis which is a consequence of the first order transition, when 
we follow the evolution of the order parameter p c coming from the superfluid phase or from the normal phase. 

3. The lower panel, corresponding to a larger hopping J = 0.3, shows a direct first order phase transition between 
the superfluid and the glass phases. Indeed, in this case the free energy of the superfluid crosses directly the 
free energy of the glass at p c . At the transition the condensate fraction jumps to zero while qEA jumps to a 
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FIG. 5: Phase diagram for the quantum model in the (fi/V, J/V) plane, for c = 3 and I = 1, for three different temperatures 
PV = 8, 15 and 30. The lines divide the phase diagram into three main regions where the system is found in a glass, superfluid 
or (normal) liquid phase. Solid lines and circles indicate the second order superfluid transition. The large dots mark the 
tricritical point where the superfluid line changes from second to first order. Dash-dotted lines and squares represent the first 
order transition between the superfluid phase and the normal (glass or liquid) phase. Dotted lines and triangles indicate the 
Kauzmann glass transition Jk(aO- The dot-dot-dashed black line indicates the normal liquid-superfluid transition for hard core 
bosons at zero temperature. Since below this line p = at T = 0, the interaction is not relevant around it, so this line is the 
T — limit of the normal liquid-superlfuid transition also for the model investigated here. The dashed black line serves as 
a guide to the eye: it has been obtained by interpolating the large fj, superfluid lines at different temperatures (see text for 
details). 



non-zero value. Also in this case we observe a phenomenon of hysteresis, and the exact first order transition 
point has to be determined by the crossing of the RS free energy of the superfluid phase and the 1RSB free 
energy of the glass phase. 

In the next section, we investigate the evolution of these phase boundaries by varying hopping, temperature and 
chemical potential. 



C. Phase diagram of the quantum model 



The quantum cavity method (at least in its present formulation) allows only to access the finite temperature 
properties of the system. We therefore start the discussion of the phase diagram from finite temperature, and then 
we discuss how to extrapolate the results to the T = limit. 



1. Finite temperature phase diagram 

The phase diagram, at finite temperature in the (//, J) plane is shown in Fig. 5. For any fixed (and low enough) 
temperature, there is a curve /z c (J), or J c (/i)i which separates a superfluid phase from a normal phase. A tricritical 
point divides this transition line into two parts. On one side, at small chemical potential, the transition is second 
order. Beyond the tricritical point instead, the transition becomes first order. This happens at bigger values of /i, 
close to the glass transition. We chose different colors to represent the phase diagram at different temperatures. 
Solid lines and circles indicate second order superfluid transitions while their continuations, dash-dotted lines and 
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FIG. 6: Left Panel: Schematic (J/V, p) phase diagram of the model. The full line represents the second order superfluid 
transition J c (p), separated by a tricritical point (large dot) from a first order superfluid transition accompanied by phase 
coexistence (dot-dashed lines). The dotted line represents the Kauzmann transition Jk(p)- In the limit T — >• 0, the transition 
lines have distinct behaviors: the first order transition line has J ^> T, and it has a finite limit for T — > 0. On the contrary, 
the other lines have J oc T and they shrink to the J = axis for T — > 0. Therefore, at T — and for J > 0, the low-p part of 
the phase diagram contains the superfluid phase while the large-p part contains the glass. The red lines indicate the behavior 
of the superfluid-glass transition at T = 0. Right Panel: Data for the (/3J, p) phase diagram, for c = 3 and I = 1, and for three 
different temperatures fiV = 8, 15 and 30. The three regions, glass, superfluid and normal liquid phase, are here reported. The 
transition lines are plotted using the same styles as in Fig. 5. As a consequence of the first order phase transition a region of 
phase coexistence - delimited by dash-dotted lines - is present. The inset shows the re-entrant behavior of the glass transition 
line, at low enough temperature. We plot ft J on the vertical axis in order to show that the transition lines are proportional to 
T in the limit T -> 0. 

squares, represent first order transition between a superfluid phase and a normal liquid or glass phase. The liquid-glass 
Kauzmann transition, pk{J) or Jr-(/x), is reported using a dotted line and triangles. It hits the horizontal axis, at 
J = 0, in correspondence of the classical glass transition; on the other side, it crosses the first order supcrfhiid-normal 
phase transition. The dynamic transition line is very close to the Kauzmann transition so we don't report it for 
clarity. 

We now describe the phase diagram as a function of the density. Since the cavity method works necessarily in the 
grand-canonical ensemble, we are forced to measure the density at fixed pi. Therefore, the density has some fluctuations 
due to numerical noise. Unfortunately, for this model the interesting part of the phase diagram is contained in a very 
small interval of density, and for this reason the noise is important and prevents us to obtain transition lines as clean 
as in the (J, p) phase diagram of Fig. 5. For this reason, in Fig. 6 we report a schematic (J, p) phase diagram together 
with the actual data. We keep the same code of colors and lines as in Fig. 5. 

The main difference between the (J, p) and (J, p) phase diagrams is that in the latter case, as already discussed, we 
observe a phase coexistence between the superfluid phase and a normal (liquid or glass) phase. The region of superfluid- 
glass phase coexistence is particularly interesting since it represents a region which manifests both amorphous order 
and off-diagonal long range order, even if they are phase separated. Besides the superfluid transition in Fig. 6 we 
reported the glass (Kauzmann) transition (dotted line) which is also expanded in the inset. From the inset it is clear 
that the glass transition as a function of the "quantum parameter" f3J, has re-entrant behavior. In fact, looking at 
the /3 = 30 or /3 = 15 curves, one can note that the system reaches the glass phase at lower densities when we switch 
on quantum fluctuations. This phenomenon, which at first sight can appear very surprising, has been recently found 
in a related work, focusing on the description of the quantum glass transition from the point of view of a microscopic 
theory of the system dynamics [35] . Moreover it has also been seen in another recent work, in the context of a simple 
optimization problem [34]. In general, one can imagine that, beyond the details which pertain to this model, this 
re-entrant behavior may have a more general interpretation, in terms of an order-by-disorder mechanism which is 
induced by quantum fluctuations and in which a particular order is selected for entropic reasons, as we will argue in 
the following. 
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2. Zero temperature limit 

We now give some hand- waving arguments on the behavior of the transition lines for T — > 0. Although we don't 
expect our arguments to be completely convincing, the following is the only consistent scenario we were able to 
elaborate that is in agreement with the observed behavior of the lines at low temperatures. 

The evolution of the phase transition lines in the {fi, J) plane of Fig. 5, suggests that in the limit T — > both 
the second order superfluid transition line and the Kauzmann transition line shrink to the origin of the axis. Note 
that we are not interested in the region \x < since in this region, at T = 0, the system has a small density and it 
behaves very similarly to weakly interacting hard core bosons. Note also that in the classical case J = 0, we know that 
is proportional to T at low temperature, which is consistent with the shrinking of the line pk{J) to the origin. 
Therefore, for T — > the glass phase "invades" the lower part of the phase diagram in Fig. 5, while the superfluid 
phase "invades" the upper part, and the first order superfluid-glass transition line at p > extends down to the origin 
(J = and p = 0). 

Unfortunately, our 1RSB code becomes very slow and unstable for too large (3 J, preventing us from drawing the 
superfluid-glass line for J » T. Still, we observe that the lines at higher temperature seem to be close to the 
continuation of the lines at lower temperature. This indicates that, as expected, for large enough J and p, the system 
has reached its zero temperature limit. At the same time, this allows us to "extrapolate" the first order superfluid- 
glass transition line at T — by taking, at each T, the largest values of J that we can access, and interpolating these 
values. The result is shown as a black dashed line in Fig. 5, and we believe that the extrapolation is a very reliable 
representation of the T = line. Therefore, at strictly zero temperature (T = 0) and positive hopping (J > 0), the 
phase diagram contains only a superfluid and a glass phase separated by a first order transition. On the other hand, 
at strictly zero hopping J = and zero temperature, one recovers the classical model of [21], which displays a RFOT 
between a liquid and a glass (however, one has to rescale p by the temperature in the classical limit, therefore the 
transition happens at p = in Fig. 5). Hence, at T = 0, the limit of vanishing hopping is extremely singular, the 
behavior of the quantum model at any J > being completely different from that of the classical model at J = 0. 

It is interesting to understand in more details how this singularity develops in the limit T — > 0. Indeed, at large (3, 
the normal liquid phase exists only in a region < p < 1/(3 and < J < 1/(3, see Fig. 5. Since both J and p go to 
zero with temperature, in this part of the phase diagram V is much larger than any other energy scale and one can 
consider it as infinite. Then, only the three energy scales J, p, T remain and the phase diagram must be a function of 
(3J and (3p only. In the limit (3 — > oo, also this energy scale disappears, and we conclude that at T — 0, J c {p) oc p, 
i.e. the superfluid-glass transition line must be linear at small J and p. 

When we eliminate the chemical potential and look to the plane (p,J), the scenario described above leads to a 
phase diagram characterized by two distinct regimes, as shown schematically in the left panel of Fig. 6: 

• If J remains finite while T — > (i.e. if (3J 3> 1), the normal liquid phase disappears. Only the superfluid 
and glass phase survive, and they are separated by a phase coexistence region, determined by the Maxwell 
construction, and delimited by lines that reach a finite limit when T — > 0. We expect that when J/V > 1 
the system is superfluid at all densities while when J/V < 1 the glass phase appears at large enough densities: 
therefore the coexistence line must have the shape reported in the schematic plot (left panel) of Fig. 6. Although 
our 1RSB code is unable to access the region J ~ V >• T, we could check that these lines have the expected 
behavior at least within the RS approximation, which should be a good approximation at least for not too small 
temperatures (we don't report these data to avoid confusion in the figures). 

• On the contrary, in the interval J < 1/(3 the system is still sensible to finite temperature effects, but as we 
argued above, we expect these effects to depend only on the quantity (3 J at fixed p, since no other energy 
scale is relevant here ( J <~ T <C V, therefore V is infinite and disappears from the problem). In this regime a 
normal liquid phase, a second order superfluid transition J c {p) and a glass transition Jk(p) still survive, but 
both critical values of J vanish proportionally to the temperature in the limit T — > 0, therefore confining the 
liquid phase to a smaller and smaller region of the phase diagram which at T = reduces to the classical region 
J = alone. 

• The phase coexistence boundaries start from the tricritical point, located at J oc T for T — > 0, and they extend 
into the large J region of the phase diagram where they have a finite limit for T — > 0. The only possibility to 
match the two regimes is that the lines become vertical in the (p, J) plane in the region (3 J <~ 1; this should be 
evident from the schematic plot in the left panel of Fig. 6. In this way we define the two values of density that 
delimit the coexistence region at T = and small J > 0. 

Looking at the data at finite temperature in right panel of Fig. 6 one sees that the superfluid transition far from 
the glass satisfies well the scaling with (3 J for p < 0.3, where the transition remains second order for all values of (3. 
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FIG. 7: A schematic plot of the typical behavior of the complexity as a function of the internal entropy of the states for various 
densities p increasing from top to bottom. 

However Fig. 6 shows that around densities of the order p ~ 0.3 the transition lines do not scale very well. To observe 
the scaling with (3 J one has to go to lower temperatures which are not easily accessed with our method. Still, the 
phase coexistence lines become almost vertical at large /3J around p <~ 0.5 which is an estimate of the coexistence 
density for T = and J > 0, according to the argument above. The glass transition line is strongly re-entrant at 
the lowest temperature (we recall that for the classical model px — 0.5725), which is consistent with this estimate. 
This means that in presence of quantum fluctuations a glass can be formed at zero temperature for densities as low 
as p = 0.5, where the corresponding classical system is in the liquid phase and quite far from the glass. 

3. Argument for the re-entrance of the glass transition line 

The re-entrance of the glass transition shown in the inset of Fig. 6 is rather unexpected and it deserves a more 
detailed discussion. It shows that if we consider the classical model at a density slightly below the Kauzmann 
transition, then it is possible to make it condense into a glass phase just by switching on quantum fluctuations, for 
low enough temperatures. This means that the quantum dynamics induced by the hopping between neighboring sites 
has a profound effect on the thermodynamical properties of the systems. In fact, the kinetic exchange selects as 
equilibrium states some states which are no more exponentially numerous, inducing the vanishing of the complexity 
function. 

One can devise an explanation of this phenomenon in terms of the entropy of the states [34] . A similar explanation 
which does not explicitly make use of entropy has been proposed in [35] . Since the re-entrance happens close to the 
classical glass transition, we can for this discussion consider V ^> fi, J, T, i.e. V = oo going back to the original 
hard sphere model of Biroli and Mezard [21]. Let us focus first on the classical model at low enough temperatures. 
Since V = oo, only zero energy configurations contribute to the measure. Our argument starts from the consideration 
that beyond the dynamical transition pd, the states (or "clusters") in which the Gibbs measure splits, are generically 
characterized by a distribution of internal entropies. The zero-temperature complexity £(s) is the logarithm of the 
number of clusters that have entropy s. Since the entropy counts the number of configurations that belong to a given 
cluster, the clusters with larger entropy have more configurations, which means that particles are more "mobile" 
inside those clusters. A second remark is that the complexity £(s) is a decreasing (in the branch of interest) convex 
function. Fig. 7 shows its typical behavior for systems belonging to the same universality class as the Biroli-Mezard 
model [22]. While many small states appear in the measure, states with higher entropies are less numerous. The 
external parameters, such as temperature and density, fix which are the dominant (equilibrium) states and then which 
is the equilibrium entropy. 

Next we assume that when quantum fluctuations are introduced, if J is sufficiently small, the distribution S(s) is 
not strongly modified. This is due to the fact that configurations belonging to different clusters are at a distance 
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corresponding to O(N) particle movements. It is not possible to go "continuously" with single particle hoppings 
from one cluster to the other, because this requires global rearrangement of the system. If one thinks in terms of 
perturbation theory in the hopping term, this happens only at extremely high orders J oi - N K Along this line of 
thought, one can safely diagonalize the Hamiltonian restricted to each state when V = oo ^> J. The major effect 
of quantum fluctuations in this regime is rather to induce, for each cluster, a kinetic energy gain proportional to its 
entropy This gain comes from the fact that higher entropy clusters are larger, therefore in those clusters particles 
can delocalize more, allowing to lower the kinetic energy. Conversely, in small entropy clusters, particles are tightly 
packed and cannot delocalize to lower their kinetic energy. In other words, the entropy of each cluster measures the 
number of "neighboring" configurations belonging to the state, namely configurations which can be reached one from 
the other through single particle movements. 

The shift of the quantum energy, proportional to J times the classical entropy, has an effect on the selection of the 
equilibrium states: it will favor states with bigger entropy because their energy is lowered more. Since these states 
are less numerous, increasing the hopping J, one expects that the complexity will be lowered, and a condensation 
transition will be eventually induced at large enough J. All these arguments are not modified when besides the 
entropy each state is also characterized by a classical energy (i.e. for finite V). A distribution of energies of the states 
must be included, but this does not change the result, just implying a more complicated temperature dependence. 
Given the complexity of the model under consideration these arguments are simply a qualitative explanation of what 
happens. However we have recovered analytically the same results "quantizing" a simpler model whose classical phase 
diagram presents the same properties of the model considered here [34] . 



D. Quantum dynamics 

From a classical point of view the dynamics of glassy systems has received a lot of attention and their slow dynamic 
behavior has been the subject of a vast part of the literature (see e.g. [1, 2]). Since dynamics and thermodynamics 
are inevitably intertwined in quantum mechanics, the study of time dependent equilibrium correlations is particularly 
interesting to understand to what extent the phenomenology of glassy quantum systems resembles that of their 
classical counterpart. In this section we analyze the local Green function and the local density-density correlation in 
imaginary time. 

1. Green function and time scales 

The Green function is defined, for — (3/2 < r < /3/2, by 

GHt) = e{r)Gl{r)+e(-T)G\(T) = ^(^(0)) , (65) 
with the advanced and retarded Green functions 

Gl(r) = (a,(r)at(0)) = ^[e-^a^a}] = I £ e^^^al^ 

a, 6 

/ V 1 I ^ 

G^r) = (aj(0)a,(r)) = ±Tr[ e -W a\e^ ' a t ] = ± £ e -V+r)E a+ rE b m&i \^ , 

a, b 

where the many-body eigenvalues and eigenstates satisfy H\"4> a ) = E a \i[) a ). The way to compute time-dependent 
correlation functions such as the above within the cavity formalism has been detailed in [41]. 

The plot in Fig. 8 shows the advanced Green function G> (r) in the supcrfluid phase and in the normal liquid phase, 
just before the Kauzmann glass transition and slightly above the dynamical transition. The Green function follows 
three main regimes. First, for r — > + one finds a fast decay of the correlations. The energy scale which fixes the 
particle decay is the interaction potential. This is the dominant contribution to the average energy of the state a\\ip a ) 
just after the insertion of one particle on a low energy state \tp a )- On the other hand the contributions from the 
chemical potential and the kinetic energy are negligible with respect to V. Then for r — > + the process brings the 
system to highly excited states first. On the contrary, on the side of holes, when r — > the interaction V doesn't 
play any role. Removing a particle cannot lead to a cost in terms of V, and for short time the dominant energy 
excitation is measured by a loss of chemical potential and the proper time scale is fixed by \i. Finally, in between 
the two exponential relaxations, the Green function present a third regime in which the decay is slower. Here in the 
region where the system is not superfluid we expect an exponential decay, with a small exponent, of the order of the 
hopping. 
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FIG. 8: Time-dependent correlation functions obtained with the RS cavity method. Upper panel: Imaginary time advanced 
Green function G>(r), for c = 3 and t = 1, at f3V = 30 and J/V — 0.1, in the superfluid phase at fi/V = 0.275 and in the 
liquid phase just before the Kauzmann transition, at n/V = 0.388. The dotted line indicates the square of the expectation 
value of the bosonic operator (a) 2 corresponding to the superfluid phase. Note that the Green function G{t) is obtained by 
the periodic image of period j3 of this function for r 6 [— /3/2, /3/2]. Bottom panel: Time dependent density-density correlation 
for the two same regimes. The dotted line indicates the value of the qEA for the normal liquid (obtained with the 1RSB cavity 
method). 



2. Density correlation and ergodicity breaking 

The second time dependent correlation which is interesting from the point of view of glassy dynamics is the density- 
density correlation, defined in imaginary time as follows: 

(n l (r)n l (0)) = ^[e^-^ h t } . (67) 

In the limit of large time and small temperatures this quantity (in particular its connected component) decays to the 
overlap parameter qEA defined in (59). In fact, rewriting Eq.(67) we obtain 

(fti(r)ni(0)) = |^e-^-^- E «)|^ |n# 6 )| 2 

a.b 

i i (68) 

The second term in (68) goes to zero when j3 J r J 1 (assuming that J < V, fj, since this is the interesting regime 
for the purpose of this discussion) , therefore in this limit 

a 

The same result is obtained for real time correlations in the large time limit, for any /?, if one assumes that terms in 
the sum with a ^ b vanish due to fast oscillations. In the following we focus on imaginary time correlations in the 
regime j3J 3> tJ 3> 1, but the same arguments can be repeated for real time correlations at any /3 and large times. 
Moreover, in the following we implicitly assume that the thermodynamic limit ./V — > oo is taken before any other limit. 
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Let us discuss first what happens in the liquid phase. Here, we expect that the connected correlation functions 
vanish in the large time limit (since the liquid phase is assumed to be ergodic). Therefore we expect that the large 
time limit of (n,(r)ni(0)) equals (hi) 2 , which leads to 



-pE a 



\(tpa\ni\tp a }\ 2 = 



-f3E a 



{lpa\hi\^ a ) 



(70) 



Defining 



Vp[m\ 



5[(tpa\hi\lpa) 



rii 



(71) 



as the probability distribution of nf over eigenstates sampled at temperature 1/(3, Eq. (70) can hold if and only if 



-> S [m 

JV-s-oo 



(ni)} 



(72) 



Since the sum is dominated by eigenvectors that have energy E a of the order of the average energy at temperature 
/?, the equation above states that the quantity nf = (V'ol^ilV'o) does not fluctuate from eigenvalue to eigenvalue in 
this range of energy. Indeed, even if up to now we chose the canonical ensemble, this requirement is unessential and 
one could have equivalently stated this concentration property of local observables also at the level of microcanonical 
averages. This assumption is known as a weak form of the eigenstate thermalization hypothesis [54-58], see [59] for a 
more detailed discussion. 

We now repeat the same discussion in the glass phase. To this aim, we need to make an important assumption on 
the spectrum of the Hamiltonian in the 1RSB phase. Indeed, as we previously discussed, the cavity method allows 
us to establish that the imaginary-time path Gibbs measure constructed via the Suzuki- Trotter formalism undergoes 
a clustering transition as in the classical case. Yet, imaginary time paths are abstract objects, and the consequences 
of this clustering transition on the spectrum of the Hamiltonian are not clear, since a priori clustering might depend 
on the particular basis that is chosen and on other details of the Suzuki- Trotter decomposition. In the following, we 
assume that in the clustered phase, the relevant eigenstates of the Hamiltonian (those which dominate the canonical 
sum) are also arranged in disconnected clusters corresponding to the thermodynamic states [9, 60]. Other states of 
the Hilbert space do not belong to any cluster, but their energies are extensively different from the thermodynamic 
energy and therefore they are exponentially suppressed in the canonical sum. We can then decompose the canonical 
sum as a sum over the clusters and a sum over the eigenstates belonging to the same cluster. A reasonable hypothesis 
that holds in the classical case and that we believe to hold also here, is that the dynamical behavior of the states 
inside each cluster is thermal (or ergodic). In other words, we assume that the reasoning we made for the liquid can 
be applied within any cluster a, leading to 



1 



iV-^oo 



(«<>a = 



-PE a 



{-ipa\hi\^ a } 



(73) 



and 



N- 



(ni) c 



(74) 



On the other hand, we expect (hi) a to fluctuate from cluster to cluster, therefore eigenstate thermalization can be 
assumed to hold at most inside each cluster, signaling, globally, a breakdown of ergodicity. Substituting (73) in (68) 
one obtains: 



(75) 



This reduces to the average of the square local density over the amorphous states. Averaging over different sites and 
subtracting the connected term one recovers the definition of Eq. (59), since W a = Z a /Z: 



(ni(r)ni(O)) - (thY 



0J>tJ>1 Z 



-i 2 



= Qea 



(76) 



28 



When the quantity in (76) goes to zero, it means that (72) holds in the whole Hilbert space and thus the system 
is ergodic. We note that a consequence of this discussion is that in the quantum case, the RS solution can give 
informations about clustering. Indeed, in the region between the dynamical and thermodynamical transitions, the 
Parisi parameter is m = 1 and in that case even in presence of clustering, local quantities can be computed using 
the RS solution. Therefore one can compute (ni(r)nj(O)) using the RS cavity method, and detect the presence of 
clustering using Eq. (76). This is not surprising, since already in the classical case the investigation of dynamics at 
the RS level can give information about clustering and 1RSB [18]. Similar results have been discussed in [29, 30]. 

In Fig. 8 we have plotted the time dependent density-density correlation for the same parameters used for the 
Green functions, namely ft = 30, J = 0.1 and two different densities such that the system is in the superfiuid and 
in the liquid phase close to the glass transition. The plot shows that the function within the glass does not relax 
at large times. It remains at a high plateau, higher than qEA- The discrepancy is due to the fact that ft J is not 
sufficiently high to allow the complete relaxation of the correlations. This is seen also in the superfiuid phase where 
the correlation remains always above zero. 

E. Exact diagonalization results 

To conclude this section, we report some results of exact diagonalization that we performed on this model. The 
aim of this investigation was mainly to check that the glass phase is not superfiuid: indeed, with the cavity method 
we do not find any non-trivial 1RSB solution with (a) ^ 0, but we cannot reach very low temperatures so we cannot 
exclude that a tiny condensate fraction might appear at a very low temperature. Moreover, with the cavity method 
we cannot access the region of the glass phase at ft J > 1 where superfluidity might be expected, because the 1RSB 
code becomes extremely slow and unstable in that region. Exact diagonalization can give in principle access to the 
exact ground state of the system, but the sizes one can study are very small. 

There is however one very important problem: in the glass phase, we expect that the ground state is highly 
degenerate in the thermodynamic limit, but at finite N, this degeneracy should be lifted leaving a unique ground 
state and a band of very low energy almost degenerate eigenstates. Defining, on a given small graph, the glass ground 
states is a very difficult task since they are well defined only in the thermodynamic limit, while at finite N they are 
linear combinations of the ground and low-energy eigenstates. Let us make our expectations more precise. Based on 
the arguments of [9, 34, 60], and on the discussion of section IV C 3, we expect that in the glass phase, at very large 
(but finite) N, one can construct almost exact low-energy eigenstates \GS a ) of the quantum Hamiltonian, each of 
them being a linear combination of classical configurations belonging only to one classical glass state (or "cluster"). 
Indeed, since different glassy states have Hamming distance of order N in the Hilbert space, we expect tunnelling 
between them to be of order ~ e"^. However, exactly because of this tiny but finite tunnelling rate, the exact 
eigenstates of the Hamiltonian at finite ./V are linear combination of the states \GS a ), and are spread in an interval of 
energy <~ e~ N above the exact ground state [60]. The lowest of them might well be an almost uniform superposition 
of all the glassy states \GS a ) (think to the classical example of a double well potential with a finite barrier), and 
in that case the local density would be uniform in the ground state making the detection of glassy order impossible 
without looking to the whole band of almost degenerate states. 

When N is very small, however, the situation might be very different. Indeed, in this case we expect that the 
number of glassy states will be quite small, since at high density on a very small graph there will be not so many 
ways of packing particles. If a given graph has too many zero-energy configurations, then at such small N particles 
will be extremely mobile and delocalizcd, and the system will be physically closer to the superfiuid than to the glass 
phase. Conversely, if there are only a few good packings, they will be arranged in clusters that are not completely 
equivalent: one of them will have the largest classical entropy (compare the discussion of section IV C 3) and the 
corresponding state \GS a ) will have lower energy than all others \GS a '). The probability of having two equivalent 
"clusters" that could mix to form an almost uniform ground state should be quite small for small values of N. To 
check this hypothesis we measured the parameter q av (recall that we can't access Qea) defined in Eq. (60), in the 
ground state for some values of the parameters deep in the glass phase (we give more details below). We found that 
g av has a strongly bimodal distribution over the random graphs: for some graphs g av is extremely small, while for 
some others it is of the order of the qEA computed by the cavity method. We found that the fraction of the latter 
graphs increases with N in the regime of small N we can access. We interpret this as a confirmation of the discussion 
above: although for N — > oo one must find q av = with probability one, the physical mechanism responsible for this 
is the huge degeneracy of equivalent clusters, which is still not present for such small N. Based on this, we assume 
that the ground state of those graphs that have a large q av is concentrated on a single glass state, and we focus on 
them for our analysis, while we disregard those graphs that have q av ~ 0. 

We stress that the above procedure is quite arbitrary and strongly biased by our expectations and by the comparison 
with the cavity method. Therefore, the results of exact diagonalization might be misleading and shall be taken with 
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FIG. 9: Ground state density (a) and g av (b) as functions of fi/V, with c = 3 and £ = 1, at hopping J/V = 0.3, obtained by 
exact diagonalization in the grand-canonical ensemble. The maroon triangles show the result of the cavity method for the same 
parameters at j3V = 15. In correspondence of the jump there is the first order superfluid-glass phase transition. Upper panels: 
The black circles, blue squares and red diamonds represent the ground state density and g av for some selected realizations of 
random graphs of size, respectively, N = 20, 22, 24. Lower panels: Average of the same quantities over 50 random graphs. 



extreme care. On the other hand, the cavity method gives the exact solution of the problem in the thermodynamic 
limit (at finite temperature), so it is reasonable to use it as a guide to interpret the exact diagonalization results. 
Indeed, a systematic study of this model only based of exact diagonalization is impossible given the technical difficulties 
discussed above and the complexity of the spectrum of the Hamiltonian (see also [34]). In the following, we give the 
details of our procedure and the results, so that the reader will be able to build her/his own interpretation of the 
data. 



1. Results in the grand-canonical ensemble 



Another technical problem in doing exact diagonalization on this system is that the Hamiltonian conserves the 
number of particles. Therefore, its eigenstates are such that (a) = 0, and moreover, the density, which is the control 
parameter, has jumps of order 1/JV at some values of the chemical potential. The latter is a particularly severe 
problem, because all the relevant transitions happen in a very small range of densities, of the order of 1/N for the 
values of N we could access. A way out of both problems is to perform exact diagonalization in the grand-canonical 
ensemble, introducing a small field h coupled to the condensate, by adding to the Hamiltonian a term — /ty^foj + a\). 
In the superfluid phase the symmetry is spontaneously broken, therefore if one takes the limit h — > after N — > oo, 
the average of (a) is finite, as we discussed in Eq. (61). In practice, one can scale the field h as h = ho/VN upon 
increasing N, which would give the correct result in the thermodynamic limit. 

We performed exact diagonalizations in the grand-canonical ensemble at fixed J = 0.3, with h — 0.09. We first 
fixed [i = 1.6 deep in the glass phase and investigated the distribution of g av over the choice of the random graph. 
We always found a bimodal distribution with a dip around q av =0.1. Therefore we discarded the graphs that had 
a value of q av < 0.1 at this fi. Table I shows the statistics of the discarded graphs that we found analyzing different 
sizes N = 20,22,24. 

We then performed, for the accepted graphs, a scan in fi to study the glass transition. In Fig. 9 we show a comparison 
between the ground state density and q av obtained with exact diagonalization and the same quantities obtained with 
the cavity method at /3 = 15. Due to the first order phase transition in /i, both have a jump at a value /i c . The upper 
panels of Fig. 9 show the density and q av for some representative graphs, while the lower panels show the average 
over 50 accepted graphs. The agreement between the exact diagonalization and the cavity method is very good in 
the RS region, while it is poor within the glass phase. The reason is that here, on the one hand exact diagonalization 
results suffer of strong finite size effects, and on the other hand the predictions of the cavity method are computed 
at finite temperature. The transition predicted by the cavity method at j3 — 15 appears at /i c ~ 1.23, but a slight 
residual dependence of fj, c on (3 is present, and if we extrapolate the results of the cavity method computed at (3 = 8 
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N 


Total graphs 


Discarded 


Accepted fraction 


20 


1470 


1420 


0.034 


22 


506 


456 


0.099 


24 


206 


156 


0.243 


26 


89 


39 


0.562 


28 


77 


27 


0.649 



TABLE I: Statistics on the number of discarded graphs (having g a v < 0.1 in the ground state) at different N and J/V = 0.3. 
For N = 20, 22, 24, we did grand-canonical computations with p/V = 1.6. For N = 26, 28, we did canonical computations with 
15, 17 occupied sites (corresponding to p — 0.577, 0.607 respectively). 




FIG. 10: Exact diagonalization results in the canonical ensemble for c = 3 and I = 1, for three different sizes N = 20, 24, 28, 
with N p = 12, 14, 17 occupied sites respectively (p ~ 0.6), as functions of J/V. Results are shown for some representative 
graphs chosen according to the criterion discussed in the text, (a) The condensate density is seen to decrease with N at small 
J, while it is almost independent of N at large J. In the inset we show the average (over 50 accepted graphs) of the largest 
eigenvalue n max of the one-body density matrix for two different hopping strenghts in the superfiuid and glass phase, as a 
function of N. The former increases linearly with N, the latter is independent of N. (b) Behavior of g a v for the same graphs. 



and (3 — 15, for J — 0.3, we obtain fi c ~ 1.3 at T = 0. The results of the exact diagonalization instead support a 
transition at a slightly larger value of p. However, the extrapolation of the results obtained by exact diagonalization 
within the grand-canonical ensemble can be non-monotonic in the size of the system especially in the glass phase 
where the average density at fixed p does not coincide for graphs of different sizes and is basically fixed by the allowed 
filling fractions. The calculation of the condensate fraction through Eq. (61) turned out to be unreliable even at small 
densities for such small sizes, so we do not report it. 



2. Results in the canonical ensemble 



An alternative way of avoiding the problem of varying the density is to work at fixed density and look at the 
superfluid-glass transition as a function of the hopping J, in the canonical ensemble. The advantage is that we can 
access slightly larger sizes in this case. We repeated the same procedure to select the graphs. For TV = 20, 22, 24, we 
kept the same graphs of the grand-canonical study. For larger N, we fixed the hopping at J = 0.3, and the density 
as close as possible to p = 0.6, and we again discarded the graphs that had a value of q av < 0.1. Tabic I shows the 
statistics of the discarded graphs for different N. 

The condensate density p c can be conveniently estimated in this case by the maximum eigenvalue n max (divided 
by N) of the one-body density matrix &ij — (ajfij), where the average is taken over the ground state [53]. Indeed, in 
presence of Bose-Einstein condensation one has p c = limTv-j-oo n max /N, while all the other eigenvalues are finite for 
N — > oo. This is a well-defined criterion for Bose-Einstein condensation in interacting systems, proposed in [53]. The 
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reduced density matrix cry is an Hermitian single particle operator that can be diagonalizcd: 

* = 5> A |A)<A| , (77) 

A 

and its spectrum {n\} can be interpreted as the distribution of the average occupation number of the single particle 
states | A). From the properties of semi-positivity and of the trace Tr[<r] = N p , where N p is the number of particles in 
the system, each eigenvalue is subject to the constraint < n\ < N p . The manifestation of Bose-Einstein condensation 
appears as an extensive occupation of a single particle state, namely n max — 0(N) and it gives an estimate of the 
condensate density of particles in that state as p c — n max /N. The presence of an extensive eigenvalue is related to 
the appearance of off-diagonal long range order. In fact, considering the limit 

lim (ataj-) = (a 4 >(%> -> n max \^ ax Xf ax , (78) 
K— i|— *-oo 

where A™ a:E = (i\\ max ), one recovers the relation between the expectation value of (hi) and the condensate density: 

^Ei^)i 2 = VEK^)i 2 = V = ^- < re > 

i i 

In Fig. 10 we show the behavior of the c/ av and of the condensate density for TV = 20, 24, 28, as a function of J at 
fixed density. The comparison of the results for the three sizes that we show in Fig. 10 is not straightforward because 
they correspond to slightly different densities. Given the small sizes that are tractable by exact diagonalization it is 
not possible to change smoothly the density. However one can note that on average n max grows linearly with TV for 
large J, while it is independent of N inside the glass phase, see the inset of Fig. 10(a). We therefore conclude that 
the exact diagonalization data support the result that the glass phase is non-superfluid even at zero temperature and 
for large hopping. 



V. CONCLUSIONS 



In this paper we studied the quantum version of a prototypical mean field lattice model of glass forming liquid, 
the Biroli-Mezard model [21]. With respect to previous studies of quantum versions of mean field glass models [27- 
33], this model has several important characteristics in the classical case: i) it presents a glass transition at zero 
temperature as a function of density; and ii) it presents a large degeneracy of glassy states with different entropies. 
We previously reported on the study of a model with these properties [34] , which, although very instructive, is quite 
abstract. The Biroli-Mezard model is instead a realistic model of interacting bosons, which can be quantized adding 
an hopping term and therefore allowing for a supcrfluid phase. We also stress once again that the classical physics 
of the Biroli-Mezard model is that of a Random First Order Transition (RFOT), therefore it should be in the same 
universality class of structural glasses, and quite different from a standard second order spin glass transition such 
as the one of the Shcrrington-Kirkpatrick model [36] . In the latter case, the addition of an hopping term gives rise 
generically to a superglass phase [37, 38]; we investigated in this paper whether this remains true for a model in the 
RFOT class. 

We investigated the phase diagram of the model by means of the 1-step replica symmetry breaking quantum cavity 
method, that we conjecture to give the exact solution of the thermodynamics of the model. We showed that the phase 
diagram is quite rich, displaying superfluid, normal liquid, and glass phases, separated by different phase transitions. 
In particular, we obtained two unexpected and quite interesting results: first of all, we showed that at low enough 
temperature the glass transition line is re-entrant as a function of quantum fluctuations, implying that one can form 
a glass by increasing quantum fluctuations at fixed density. Similar results have been obtained in [35] by an extension 
of Mode-Coupling Theory to quantum hard spheres. Since Mode-Coupling Theory is an approximate description of 
the glass transition based on some uncontrolled approximations, it is very important to confirm the result of [35] by 
means of an exact solution of a mean field model. Additionally, we showed that the standard RFOT glass transition 
is replaced by a first order supcrfluid-glass transition at zero temperature, accompanied by phase coexistence between 
the two phases, while at the same time the glass transition completely disappears. This shows that for models with 
such a complex phase space in the classical limit, introducing quantum fluctuations has a dramatically singular effect, 
changing completely the nature of the transition between the liquid and the glass phases. Moreover, the first order 
superfluid- glass transition is accompanied by a jump in density, implying that there exists an interval of densities 
where the two phases would coexist in a finite dimensional version of the model. One would therefore obtain a 
simultaneous presence of superfluid and glassy ordering, which however would not give rise to a true superglass phase 



32 



since they would be phase separated. At variance to what happens in models displaying a second order spin glass 
transition [37, 38], we did not find any pure "superglass" phase in the quantum Biroli-Mezard model: in more technical 
terms this means that, for the range of parameters considered here, we have not found a solution of the quantum 
1-RSB equations within the broken phase (a) ^ 0. 

These predictions will hopefully bring new insights on the current theoretical debate on the nature of superfluidity 
in solid Helium [11]. Moreover, the present study is quite important, in our opinion, since it opens the way to extend 
to the quantum regime the replica techniques that have been used in the description of classical glasses [24] . These 
techniques are based on the assumption that the glass transition is a Random First Order Transition. Since the 
present study shows that this is not true in the fully quantum regime, the method of [24] will have to be adapted to 
describe the crossover between the RFOT and the first order regime. We expect that this should be possible using 
the phase diagram of the model investigated here as a guide. If successful, this program should lead to a quantitative 
computation of the glassy phase diagram for realistic potentials such as the one of Helium, which would in turn allow 
for a systematic numerical and experimental investigation, in order to link consistently the results of [8] with the 
established picture of classical glasses. 
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